PREDICCIÓN DE LOS SALARIOS TRIMESTRALES EN EL ESTADO DE NEW YORK¶

In [1]:
import pandas as pd
import numpy as np
from sklearn.preprocessing import LabelEncoder
import plotly.express as px
import matplotlib.pyplot as plt
import pmdarima as pm
import missingno as msno
import seaborn as sns
import scipy.stats as stats
from sklearn.model_selection import train_test_split
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import statsmodels.api as sm
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_percentage_error
from statsmodels.tsa.seasonal import seasonal_decompose
import xgboost as xgb
from xgboost import XGBRegressor
from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GridSearchCV

DATASETS

Enlace para descargar los datasets: https://drive.google.com/drive/folders/1hYewiVm4w3-GMIWixnP6Lc5_5HtLJYIo?usp=sharing

El objetivo de este trabajo es construir un modelo predictivo de los salarios trimestrales en el estado de New York.

Este trabajo consta de un conjunto de datos principal 'census' y cuatro conjuntos de datos abiertos. Los datos abiertos contienen información macroeconómica y geográfica y han sido agregados al trabajo con el propósito de mejorar el modelo predictivo.

CENSUS
El conjunto de datos principal de este trabajo procede del Quarterly Census of Employment and Wages program (QCEW), la principal fuente de datos del mercado laboral que existe en el estado de New York. Este conjunto de datos recopila las principales informaciones sobre empleo y salario a nivel estatal. Proporciona un censo virtual del 97% de los empleados y empleadores no agrícolas del estado, lo que lo hace una fuente de datos muy amplia y con esto precisa. Debido a su naturaleza, el conjunto de datos principal se presta a ser cumplimentado con datos abiertos relativos a variables macroeconómicas.

GDP
El dataset 'GDP' contiene los datos trimestrales (en trillones de dólares) relativos al producto interno bruto en el estado de New York de 2000 a 2022. Los valores desde el año 2000 hasta finales de 2004, han sido estimados considerando el GDP anual del estado de New York.

INFLATION
El dataset de 'Inflation' contiene las tasas de inflación trimestral en EE.UU. desde el año 2000 hasta finales de 2019.

UNEMPLOYMENT
El dataset 'Unemployment' contiene los datos de desempleo detallados por mano de obra, empleados, desempleados y tasa de desempleo. Todos estos datos están detallados por año, més y área geográfica.

ÁREA
El dataset 'Area' contiene la información relativa a las subdivisiones administrativas de las áreas geográficas en el estado de New York.

  • County. Los condados son los pilares utilizados para construir áreas geográficas cada vez más grandes para las cuales se reportan estadísticas del mercado laboral. El estado de New York cuenta con 62 condados.

  • Metropolitan Statistical Area. Cada Área Estadística Metropolitana debe tener al menos un área urbana de 50.000 o más habitantes y consta de uno o más condados asociados. El dataset cuenta con 12 Áreas Estadísticas Metropolitanas.

  • Micropolitan Statistical Area. Cada Área Estadística Micropolitana debe tener al menos un área urbana de al menos 10.000 pero menos de 50.000 habitantes y consta de uno o más condados asociados. El estado de New York cuenta con 14 Áreas Estadísticas Micropolitanas.

  • Metropolitan Division. La División Metropolitana es una agrupación más pequeña de los condados que se cumple en caso un Área Estadística Metropolitana contenga un único núcleo con una población de 2,5 millones de habitantes. El estado de New York cuenta con 2 Divisiones Metropolitanas.

  • Workforce Investment Region. Es el organismo para el desarrollo de la fuerza laboral. El estado de New York con 33 juntas de desarrollo de la fuerza laboral local. Labor Market Area. Las Regiones del Mercado Laboral (LMR) corresponden a diez áreas geográficas del estado, que están definidas por el Departamento de Trabajo del Estado de Nueva York.

In [2]:
census = pd.read_csv("C:\\Users\\giuli\\UCM\\TFM\\census_2000_2019.csv")
gdp = pd.read_csv("C:\\Users\\giuli\\UCM\\TFM\\GDP.csv")
inflation = pd.read_csv("C:\\Users\\giuli\\UCM\\TFM\\inflation_USA_2000_2019.csv")
unemployment = pd.read_csv("C:\\Users\\giuli\\UCM\\TFM\\local-area-unemployment-statistics-beginning-1976-1.csv")
Area = pd.read_csv("C:\\Users\\giuli\\UCM\\TFM\\Area.csv")

PREPROCESADO

PREPROCESADO DE CENSUS

El dataset ‘Census’ cuenta con 12 variables y 982173 observaciones.

De una primera observación se nota que la variable ‘Total Wage’ tiene valores nulos.

En primer lugar, para facilitar la gestión de los datos se divide por 1.000.000 el valor de ‘Total Wage’ y se almacena en una nueva variable ‘Wage_million$’.

Se crea una nueva variable ‘Employment’ donde se almacena la suma de la población empleada en cada trimestre.

Se necesita crear una referencia temporal trimestral de tipo datetime. Para lograrlo, se remplazan los valores 1, 2, 3, 4 de la variable ‘Quarter’ por 31-03, 30-06, 30-09, 31-12 y al final de cada fecha se agrega la información del año correspondiente, procedente de la variable ‘Year’. Al final se transforma el objeto string en objeto datetime, con su formato correspondiente.

En cuanto a la variable ‘Área’ se decide quedar con la unidad administrativa más pequeña, o sea el condado. Se eliminan todas las observaciones referentes a unidades administrativas diferentes. Por dicha razón, se excluyen de la columna ‘Área’ todas las observaciones que no contienen la palabra ‘county’.

La variable ‘NAICS Title’ deja ver que hay una agrupación de datos ‘Total, All Industries’. Esta agrupación podría distorsionar el análisis, así que se eliminan todas las observaciones a las cuales les corresponde este valor.

En fin, se eliminan las columnas superfluas ‘Total Wage’, ‘index’, ‘Month 1 Employment’, ‘Month 2 Employment’, ‘Month 3 Employment’, ‘Area Type’.

In [3]:
census.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 982173 entries, 0 to 982172
Data columns (total 12 columns):
 #   Column                                                                                                                             Non-Null Count   Dtype  
---  ------                                                                                                                             --------------   -----  
 0   index                                                                                                                              982173 non-null  int64  
 1   Area Type                                                                                                                          982173 non-null  object 
 2   Area                                                                                                                               982173 non-null  object 
 3   Year                                                                                                                               982173 non-null  int64  
 4   Quarter                                                                                                                            982173 non-null  int64  
 5   NAICS                                                                                                                              982173 non-null  int64  
 6   NAICS Title                                                                                                                        982173 non-null  object 
 7   Establishments                                                                                                                     982173 non-null  int64  
 8   Month 1 Employment                                                                                                                 982173 non-null  int64  
 9   Month 2 Employment                                                                                                                 982173 non-null  int64  
 10  Month 3 Employment                                                                                                                 982173 non-null  int64  
 11  Total Wage                                                                                                                         982145 non-null  float64
dtypes: float64(1), int64(8), object(3)
memory usage: 89.9+ MB
In [4]:
census.head()
Out[4]:
index Area Type Area Year Quarter NAICS NAICS Title Establishments Month 1 Employment Month 2 Employment Month 3 Employment Total Wage
0 0 State New York State 2019 1 0 Total, All Industries 640172 9341688 9392837 9446466 2.001011e+11
1 1 County Albany County 2019 1 0 Total, All Industries 10177 231193 231993 232076 3.439674e+09
2 2 County Allegany County 2019 1 0 Total, All Industries 912 12748 13230 13254 1.263470e+08
3 3 County Bronx County 2019 1 0 Total, All Industries 18622 317542 319956 322237 4.528452e+09
4 4 County Broome County 2019 1 0 Total, All Industries 4381 83922 85202 85620 9.918101e+08
In [5]:
census.columns = census.columns.str.strip()
In [6]:
census["Wage_million$"] = census["Total Wage"]/1000000
In [7]:
census['Employment'] = census['Month 1 Employment']+census['Month 2 Employment']+census['Month 3 Employment']
In [8]:
census["Quarter"] = census["Quarter"].astype('string')
census["Year"] = census["Year"].astype('string')
In [9]:
census['Quarter'] = census['Quarter'].replace(['1'], ['31-03-'])
census['Quarter'] = census['Quarter'].replace(['2'], ['30-06-'])
census['Quarter'] = census['Quarter'].replace(['3'], ['30-09-'])
census['Quarter'] = census['Quarter'].replace(['4'], ['31-12-'])
In [10]:
census['Quarter'] = census['Quarter']+(census['Year'])
In [11]:
census['Quarter'] = pd.to_datetime(census['Quarter'], format='%d-%m-%Y',errors='raise')
In [12]:
census= census.drop(census[census['Area Type'] != 'County'].index)
census['County']=census['Area']
In [13]:
census = census.loc[~census['NAICS Title'].str.contains('Total', case=False)]
In [14]:
census = census.drop(['Total Wage'], axis=1)
census = census.drop(['index'], axis=1)
census = census.drop(["Month 1 Employment"], axis=1)
census = census.drop(["Month 2 Employment"], axis=1)
census = census.drop(["Month 3 Employment"], axis=1)
census=census.drop(['Area Type'], axis=1)
census=census.drop(['NAICS'], axis=1)
In [15]:
census.head()
Out[15]:
Area Year Quarter NAICS Title Establishments Wage_million$ Employment County
142 Broome County 2018 2018-03-31 Florists 6 0.184196 121 Broome County
240 Essex County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 30 0.667383 291 Essex County
269 Albany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 1.971956 591 Albany County
270 Allegany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 0.930796 358 Allegany County
271 Bronx County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 4 0.213331 65 Bronx County

PREPROCESADO DE GDP

El dataset ‘gdp’ contiene 2 variables y 92 observaciones.
Para obtener la referencia temporal trimestral que sucesivamente se necesitará para incorporar este dataset al conjunto de datos principal, se transforma la columna ‘Year’ dejando solo los 11 últimos dígitos y se pasa a formato datetime.

In [16]:
gdp.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 92 entries, 0 to 91
Data columns (total 2 columns):
 #   Column         Non-Null Count  Dtype 
---  ------         --------------  ----- 
 0   Years          92 non-null     object
 1   GDP_trillion$  92 non-null     int64 
dtypes: int64(1), object(1)
memory usage: 1.6+ KB
In [17]:
gdp.head()
Out[17]:
Years GDP_trillion$
0 2022-10-01/2022-12-31 20984356
1 2022-07-01/2022-09-30 20705577
2 2022-04-01/2022-06-30 20349123
3 2022-01-01/2022-03-31 20088134
4 2021-10-01/2021-12-31 19866796
In [18]:
gdp['Quarter'] = gdp['Years'].str.slice(11)
In [19]:
gdp = gdp.drop(['Years'], axis=1)
In [20]:
gdp['Quarter'] = pd.to_datetime(gdp['Quarter'], format='%Y-%m-%d')
In [21]:
gdp.head()
Out[21]:
GDP_trillion$ Quarter
0 20984356 2022-12-31
1 20705577 2022-09-30
2 20349123 2022-06-30
3 20088134 2022-03-31
4 19866796 2021-12-31

PREPROCESADO DE INFLATION

El dataset ‘Inflation’ está compuesto por 3 variables y 80 observaciones.
Se necesita crear una columna con la referencia temporal trimestral para luego poder agregar este dataset al conjunto de datos principal. Por dicha razón, se reemplazan los valores 1, 2, 3, 4 correspondientes a los trimestres por 31-03, 30-06, 30-09, 31-12 y se le añade el año correspondiente. Todo se convierte a datetime especificando el formato correspondiente.

In [22]:
inflation.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 80 entries, 0 to 79
Data columns (total 3 columns):
 #   Column       Non-Null Count  Dtype  
---  ------       --------------  -----  
 0   Year         80 non-null     int64  
 1   Quarter      80 non-null     int64  
 2   Inflation %  80 non-null     float64
dtypes: float64(1), int64(2)
memory usage: 2.0 KB
In [23]:
inflation.head()
Out[23]:
Year Quarter Inflation %
0 2000 1 3.4
1 2000 2 1.8
2 2000 3 2.5
3 2000 4 2.2
4 2001 1 2.7
In [24]:
inflation["Quarter"] = inflation["Quarter"].astype('string')
inflation["Year"] = inflation["Year"].astype('string')
In [25]:
inflation['Quarter'] = inflation['Quarter'].replace({'1': '31-03-', '2': '30-06-', '3': '30-09-', '4': '31-12-'})
In [26]:
inflation['Quarter']=inflation['Quarter']+inflation['Year']
inflation = inflation.drop(['Year'], axis=1)
In [27]:
inflation['Quarter'] = pd.to_datetime(inflation['Quarter'], format='%d-%m-%Y')
In [28]:
inflation.head()
Out[28]:
Quarter Inflation %
0 2000-03-31 3.4
1 2000-06-30 1.8
2 2000-09-30 2.5
3 2000-12-31 2.2
4 2001-03-31 2.7

PREPROCESADO DE UNEPLOYEMENT

El dataset ‘Unemployment’ contiene 7 variables con 75712 observaciones cada una.

Para crear la referencia temporal trimestral, que luego permitirá incorporar este dataset al conjunto de datos principal, se aportan los mismos cambios ya ejecutados en los otros datasets. La variable ‘Unemployment Rate’ es mensual y debido a la necesidad de agrupar los datos trimestralmente para que el dataset sea compatible con el resto de datos, es necesario agruparla por trimestre y área, considerando los valores medios.

La variable ‘Área’ contiene algunas observaciones que se refieren a unidades administrativas diferentes de los condados, así que hay que filtrar estos valores para excluirlos del análisis y preparar el dataset para que sea compatible con el conjunto principal.

Se eliminan las variables sobrantes 'Month','Year','Labor Force','Employed','Unemployed'

In [29]:
unemployment.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 75712 entries, 0 to 75711
Data columns (total 7 columns):
 #   Column                                                                                                                                                                                   Non-Null Count  Dtype  
---  ------                                                                                                                                                                                   --------------  -----  
 0   Area                                                                                                                                                                                     75712 non-null  object 
 1   Year                                                                                                                                                                                     75712 non-null  int64  
 2   Month                                                                                                                                                                                    75712 non-null  int64  
 3   Labor Force                                                                                                                                                                              75712 non-null  int64  
 4   Employed                                                                                                                                                                                 75712 non-null  int64  
 5   Unemployed                                                                                                                                                                               75712 non-null  int64  
 6   Unemployment Rate                                                                                                                                                                        75712 non-null  float64
dtypes: float64(1), int64(5), object(1)
memory usage: 4.0+ MB
In [30]:
unemployment.head()
Out[30]:
Area Year Month Labor Force Employed Unemployed Unemployment Rate
0 Albany City 2018 10 48200 46400 1800 3.8
1 Albany City 2018 9 47600 45600 2000 4.1
2 Albany City 2018 8 47300 45100 2300 4.8
3 Albany City 2018 7 47800 45400 2400 4.9
4 Albany City 2018 6 48200 45800 2400 5.0
In [31]:
unemployment["Year"] = unemployment["Year"].astype('string')
unemployment["Month"] = unemployment["Month"].astype('string')
unemployment["Area"] = unemployment["Area"].astype('category')
In [32]:
unemployment['Quarter'] = unemployment['Month'].replace({'1': '31-03', '2': '31-03', '3': '31-03','4': '30-06', '5': '30-06', '6': '30-06','7': '30-09', '8': '30-09', '9': '30-09','10': '31-12', '11': '31-12', '12': '31-12'})
In [33]:
unemployment['Quarter'] = (unemployment["Quarter"] + ("-") + unemployment["Year"])
unemployment['Quarter']
Out[33]:
0        31-12-2018
1        30-09-2018
2        30-09-2018
3        30-09-2018
4        30-06-2018
            ...    
75707    30-06-1990
75708    30-06-1990
75709    31-03-1990
75710    31-03-1990
75711    31-03-1990
Name: Quarter, Length: 75712, dtype: string
In [34]:
unemployment['Quarter'] = pd.to_datetime(unemployment['Quarter'], format='%d-%m-%Y')
In [35]:
unemployment.columns = unemployment.columns.str.strip()
In [36]:
unemployment['Unemployment Rate'] = unemployment.groupby(['Quarter', 'Area'])['Unemployment Rate'].transform('mean')
In [37]:
unemployment = unemployment.loc[unemployment['Area'].str.contains('County', case=False)]
In [38]:
unemployment = unemployment.drop(['Month','Year','Labor Force','Employed','Unemployed'], axis=1)
unemployment= unemployment.drop_duplicates()
In [39]:
unemployment.head()
Out[39]:
Area Unemployment Rate Quarter
346 Albany County 3.100000 2018-12-31
347 Albany County 3.700000 2018-09-30
350 Albany County 3.833333 2018-06-30
353 Albany County 4.633333 2018-03-31
356 Albany County 4.066667 2017-12-31
In [40]:
unemployment.info()

<class 'pandas.core.frame.DataFrame'>
Index: 7192 entries, 346 to 74671
Data columns (total 3 columns):
 #   Column             Non-Null Count  Dtype         
---  ------             --------------  -----         
 0   Area               7192 non-null   category      
 1   Unemployment Rate  7192 non-null   float64       
 2   Quarter            7192 non-null   datetime64[ns]
dtypes: category(1), datetime64[ns](1), float64(1)
memory usage: 192.4 KB

PREPROCESADO DE AREA

El dataset ‘Área’ cuenta con 6 variables y 62 observaciones. Se observan valores nulos debidos a que algunas de las categorías de subdivisión administrativa se cumplimentan entre sí.

In [41]:
Area.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 62 entries, 0 to 61
Data columns (total 6 columns):
 #   Column                         Non-Null Count  Dtype 
---  ------                         --------------  ----- 
 0   County                         62 non-null     object
 1   Metropolitan Statistical Area  29 non-null     object
 2   Micropolitan Statistical Area  14 non-null     object
 3   Metropolitan Division          10 non-null     object
 4   Workforce Investment Region    62 non-null     object
 5   Labor Market Area              62 non-null     object
dtypes: object(6)
memory usage: 3.0+ KB
In [42]:
Area.head()
Out[42]:
County Metropolitan Statistical Area Micropolitan Statistical Area Metropolitan Division Workforce Investment Region Labor Market Area
0 Albany County Albany-Schenectady-Troy Metro Area NaN NaN Albany-Rensselaer-Schenectady Counties Capital Region
1 Allegany County NaN NaN NaN Allegany-Cattaraugus Western New York Region
2 Bronx County NaN NaN New York City Region New York City New York City Region
3 Broome County Binghamton Metro Area NaN NaN Broome-Tioga Southern Tier Region
4 Cattaraugus County NaN Olean, NY Micropolitan Statistical Area NaN Allegany-Cattaraugus Western New York Region
In [43]:
Area['LMA'] = Area['Labor Market Area'].astype('category')
Area['WIR'] = Area['Workforce Investment Region'].astype('category')
Area['County']=Area['County'].astype('category')
Area = Area.drop(['Labor Market Area','Workforce Investment Region','Metropolitan Statistical Area', 'Micropolitan Statistical Area','Metropolitan Division'], axis=1)
In [44]:
Area.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 62 entries, 0 to 61
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype   
---  ------  --------------  -----   
 0   County  62 non-null     category
 1   LMA     62 non-null     category
 2   WIR     62 non-null     category
dtypes: category(3)
memory usage: 4.5 KB

UNIÓN DE DATASETS

Tras preprocesar los datos de todos los datasets, se va a preparar un único conjunto de datos.
‘census_1’ es el dataset generado por la unión de ‘census’ y ‘gdp’ tomando como referencia la variable ‘Quarter’. En el paso sucesivo, se le añade a este dataset la información contenida en ‘inflation’, tomando como referencia ‘Quarter’. El resultado de esta unión será el dataset ‘census_2’, al cual se le agrega el dataset ‘unemployment’, tomando como referencia ‘Quarter’ y ‘Area’. Este nuevo dataset se denomina ‘census_3’. En fin se crea el dataset ‘census_4’ que incluye el dataset anterior más la información del dataset ‘Area’, tomando como referencia la columna ‘County’.

El nuevo conjunto de datos ‘census_4’ tiene 13 variables y 472.807 observaciones.

Al parecer, las variables ‘Wage_million$’ y ‘Unemployment Rate’ tienen valores faltantes por gestionar.
Además, con el método describe aplicado al conjunto de variables numéricas, se observa que el valor mínimo de ‘Wage_million’ es 0. Esto podría ser un missing no declarado. Tras comprobar que son solo 9 observaciones en todo el dataset a llevar valor 0 en esta variable se decide eliminarlas, ya que tienen muy baja incidencia considerando la totalidad del conjunto de datos.

Antes de gestionar los datos faltantes, se eliminan del dataset las variables geográficas que no se van a incluir en el análisis (‘County’,’WIR’), ya que se escoge la variable 'LMA'.

In [45]:
census.head()
Out[45]:
Area Year Quarter NAICS Title Establishments Wage_million$ Employment County
142 Broome County 2018 2018-03-31 Florists 6 0.184196 121 Broome County
240 Essex County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 30 0.667383 291 Essex County
269 Albany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 1.971956 591 Albany County
270 Allegany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 0.930796 358 Allegany County
271 Bronx County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 4 0.213331 65 Bronx County
In [46]:
census_1= pd.merge(census, gdp, on=['Quarter'], how='left')
census_1
Out[46]:
Area Year Quarter NAICS Title Establishments Wage_million$ Employment County GDP_trillion$
0 Broome County 2018 2018-03-31 Florists 6 0.184196 121 Broome County 16627803
1 Essex County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 30 0.667383 291 Essex County 17529646
2 Albany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 1.971956 591 Albany County 17529646
3 Allegany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 0.930796 358 Allegany County 17529646
4 Bronx County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 4 0.213331 65 Bronx County 17529646
... ... ... ... ... ... ... ... ... ...
472802 Sullivan County 2000 2000-12-31 Unclassified 21 0.191175 155 Sullivan County 8412000
472803 Tompkins County 2000 2000-12-31 Unclassified 10 0.104896 50 Tompkins County 8412000
472804 Ulster County 2000 2000-12-31 Unclassified 49 0.360702 298 Ulster County 8412000
472805 Wayne County 2000 2000-12-31 Unclassified 15 0.179872 52 Wayne County 8412000
472806 Westchester County 2000 2000-12-31 Unclassified 535 9.324520 3120 Westchester County 8412000

472807 rows × 9 columns

In [47]:
census_2= pd.merge(census_1, inflation, on=['Quarter'], how='left')
census_2
Out[47]:
Area Year Quarter NAICS Title Establishments Wage_million$ Employment County GDP_trillion$ Inflation %
0 Broome County 2018 2018-03-31 Florists 6 0.184196 121 Broome County 16627803 2.8
1 Essex County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 30 0.667383 291 Essex County 17529646 0.9
2 Albany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 1.971956 591 Albany County 17529646 0.9
3 Allegany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 0.930796 358 Allegany County 17529646 0.9
4 Bronx County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 4 0.213331 65 Bronx County 17529646 0.9
... ... ... ... ... ... ... ... ... ... ...
472802 Sullivan County 2000 2000-12-31 Unclassified 21 0.191175 155 Sullivan County 8412000 2.2
472803 Tompkins County 2000 2000-12-31 Unclassified 10 0.104896 50 Tompkins County 8412000 2.2
472804 Ulster County 2000 2000-12-31 Unclassified 49 0.360702 298 Ulster County 8412000 2.2
472805 Wayne County 2000 2000-12-31 Unclassified 15 0.179872 52 Wayne County 8412000 2.2
472806 Westchester County 2000 2000-12-31 Unclassified 535 9.324520 3120 Westchester County 8412000 2.2

472807 rows × 10 columns

In [48]:
census_3 = pd.merge(census_2, unemployment, on=['Quarter','Area'], how='left')
census_3
Out[48]:
Area Year Quarter NAICS Title Establishments Wage_million$ Employment County GDP_trillion$ Inflation % Unemployment Rate
0 Broome County 2018 2018-03-31 Florists 6 0.184196 121 Broome County 16627803 2.8 6.700000
1 Essex County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 30 0.667383 291 Essex County 17529646 0.9 NaN
2 Albany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 1.971956 591 Albany County 17529646 0.9 NaN
3 Allegany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 0.930796 358 Allegany County 17529646 0.9 NaN
4 Bronx County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 4 0.213331 65 Bronx County 17529646 0.9 NaN
... ... ... ... ... ... ... ... ... ... ... ...
472802 Sullivan County 2000 2000-12-31 Unclassified 21 0.191175 155 Sullivan County 8412000 2.2 4.233333
472803 Tompkins County 2000 2000-12-31 Unclassified 10 0.104896 50 Tompkins County 8412000 2.2 2.900000
472804 Ulster County 2000 2000-12-31 Unclassified 49 0.360702 298 Ulster County 8412000 2.2 3.266667
472805 Wayne County 2000 2000-12-31 Unclassified 15 0.179872 52 Wayne County 8412000 2.2 3.466667
472806 Westchester County 2000 2000-12-31 Unclassified 535 9.324520 3120 Westchester County 8412000 2.2 3.100000

472807 rows × 11 columns

In [49]:
Area.head()
Out[49]:
County LMA WIR
0 Albany County Capital Region Albany-Rensselaer-Schenectady Counties
1 Allegany County Western New York Region Allegany-Cattaraugus
2 Bronx County New York City Region New York City
3 Broome County Southern Tier Region Broome-Tioga
4 Cattaraugus County Western New York Region Allegany-Cattaraugus
In [50]:
census_4 = pd.merge(census_3, Area, on=['County'], how='left')
census_4
Out[50]:
Area Year Quarter NAICS Title Establishments Wage_million$ Employment County GDP_trillion$ Inflation % Unemployment Rate LMA WIR
0 Broome County 2018 2018-03-31 Florists 6 0.184196 121 Broome County 16627803 2.8 6.700000 Southern Tier Region Broome-Tioga
1 Essex County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 30 0.667383 291 Essex County 17529646 0.9 NaN North Country Region Clinton-Essex-Franklin-Hamilton
2 Albany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 1.971956 591 Albany County 17529646 0.9 NaN Capital Region Albany-Rensselaer-Schenectady Counties
3 Allegany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 0.930796 358 Allegany County 17529646 0.9 NaN Western New York Region Allegany-Cattaraugus
4 Bronx County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 4 0.213331 65 Bronx County 17529646 0.9 NaN New York City Region New York City
... ... ... ... ... ... ... ... ... ... ... ... ... ...
472802 Sullivan County 2000 2000-12-31 Unclassified 21 0.191175 155 Sullivan County 8412000 2.2 4.233333 Hudson Valley Region Sullivan
472803 Tompkins County 2000 2000-12-31 Unclassified 10 0.104896 50 Tompkins County 8412000 2.2 2.900000 Southern Tier Region Tompkins
472804 Ulster County 2000 2000-12-31 Unclassified 49 0.360702 298 Ulster County 8412000 2.2 3.266667 Hudson Valley Region Ulster
472805 Wayne County 2000 2000-12-31 Unclassified 15 0.179872 52 Wayne County 8412000 2.2 3.466667 Finger Lakes Region Ontario-Seneca-Wayne-Yates
472806 Westchester County 2000 2000-12-31 Unclassified 535 9.324520 3120 Westchester County 8412000 2.2 3.100000 Hudson Valley Region Putnam-Balance Of Westchester

472807 rows × 13 columns

In [51]:
census_4.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 472807 entries, 0 to 472806
Data columns (total 13 columns):
 #   Column             Non-Null Count   Dtype         
---  ------             --------------   -----         
 0   Area               472807 non-null  object        
 1   Year               472807 non-null  string        
 2   Quarter            472807 non-null  datetime64[ns]
 3   NAICS Title        472807 non-null  object        
 4   Establishments     472807 non-null  int64         
 5   Wage_million$      472779 non-null  float64       
 6   Employment         472807 non-null  int64         
 7   County             472807 non-null  object        
 8   GDP_trillion$      472807 non-null  int64         
 9   Inflation %        472807 non-null  float64       
 10  Unemployment Rate  447126 non-null  float64       
 11  LMA                472807 non-null  category      
 12  WIR                472807 non-null  category      
dtypes: category(2), datetime64[ns](1), float64(3), int64(3), object(3), string(1)
memory usage: 40.6+ MB
In [52]:
census_4.head()
Out[52]:
Area Year Quarter NAICS Title Establishments Wage_million$ Employment County GDP_trillion$ Inflation % Unemployment Rate LMA WIR
0 Broome County 2018 2018-03-31 Florists 6 0.184196 121 Broome County 16627803 2.8 6.7 Southern Tier Region Broome-Tioga
1 Essex County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 30 0.667383 291 Essex County 17529646 0.9 NaN North Country Region Clinton-Essex-Franklin-Hamilton
2 Albany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 1.971956 591 Albany County 17529646 0.9 NaN Capital Region Albany-Rensselaer-Schenectady Counties
3 Allegany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 0.930796 358 Allegany County 17529646 0.9 NaN Western New York Region Allegany-Cattaraugus
4 Bronx County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 4 0.213331 65 Bronx County 17529646 0.9 NaN New York City Region New York City
In [53]:
census_4.describe(include=np.number)
Out[53]:
Establishments Wage_million$ Employment GDP_trillion$ Inflation % Unemployment Rate
count 472807.000000 472779.000000 4.728070e+05 4.728070e+05 472807.000000 447126.000000
mean 189.544520 40.733262 8.077940e+03 1.296661e+07 1.799564 5.787943
std 772.741135 430.958390 3.727073e+04 2.981696e+06 1.426124 1.843885
min 1.000000 0.000000 2.000000e+00 8.412000e+06 -5.600000 2.433333
25% 9.000000 0.678741 2.600000e+02 1.040547e+07 1.200000 4.366667
50% 26.000000 2.573643 9.160000e+02 1.260281e+07 2.000000 5.300000
75% 87.000000 11.443832 3.691000e+03 1.584316e+07 2.700000 7.033333
max 20551.000000 48118.821276 1.227740e+06 1.789074e+07 4.500000 14.133333
In [54]:
(census_4['Wage_million$'] == 0).sum()
Out[54]:
9
In [55]:
census_4 = census_4[census_4['Wage_million$'] != 0]
In [56]:
census_4.head()
Out[56]:
Area Year Quarter NAICS Title Establishments Wage_million$ Employment County GDP_trillion$ Inflation % Unemployment Rate LMA WIR
0 Broome County 2018 2018-03-31 Florists 6 0.184196 121 Broome County 16627803 2.8 6.7 Southern Tier Region Broome-Tioga
1 Essex County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 30 0.667383 291 Essex County 17529646 0.9 NaN North Country Region Clinton-Essex-Franklin-Hamilton
2 Albany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 1.971956 591 Albany County 17529646 0.9 NaN Capital Region Albany-Rensselaer-Schenectady Counties
3 Allegany County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 18 0.930796 358 Allegany County 17529646 0.9 NaN Western New York Region Allegany-Cattaraugus
4 Bronx County 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting 4 0.213331 65 Bronx County 17529646 0.9 NaN New York City Region New York City
In [57]:
census_4.columns = census_4.columns.str.strip()
In [58]:
census_4 = census_4.drop('County', axis=1)
census_4 = census_4.drop('WIR', axis=1)
census_4 = census_4.drop(['Area'], axis=1)
census_4=census_4.reset_index()
In [59]:
census_4["NAICS Title"] = census_4["NAICS Title"].astype('category')
In [60]:
census_4.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 472798 entries, 0 to 472797
Data columns (total 11 columns):
 #   Column             Non-Null Count   Dtype         
---  ------             --------------   -----         
 0   index              472798 non-null  int64         
 1   Year               472798 non-null  string        
 2   Quarter            472798 non-null  datetime64[ns]
 3   NAICS Title        472798 non-null  category      
 4   Establishments     472798 non-null  int64         
 5   Wage_million$      472770 non-null  float64       
 6   Employment         472798 non-null  int64         
 7   GDP_trillion$      472798 non-null  int64         
 8   Inflation %        472798 non-null  float64       
 9   Unemployment Rate  447117 non-null  float64       
 10  LMA                472798 non-null  category      
dtypes: category(2), datetime64[ns](1), float64(3), int64(4), string(1)
memory usage: 33.8 MB

ANÁLISIS Y GESTIÓN DE LOS NULOS

La variable ‘NAICS Title’ presenta una categoria de ‘Unclassified’, este valor se pasa a nulo.

Se observan cuántos valores nulos hay en cada variable y su proporción con respecto al total de los datos.
En general hay pocos valores nulos. Los nulos de la variable ‘Wage_million’ son menos del 0,001%, los de ‘NAICS Title’ son aproximadamente un 1% y los de la variable ‘Unemployment Rate’ son casi un 6% del total de los datos.

Los valores nulos no tienen correlación entre sí.

Se eliminan todas las observaciones que contienen valores nulos en las variables ‘Wage_million$’ y ‘NAICS Title’ debido a su baja incidencia con respecto al total de los datos. Mientras que para la variable ‘Unemployment Rate’, los valores nulos se estiman con la media de los otros valores.

In [61]:
census_4['NAICS Title'] = census_4['NAICS Title'].replace('Unclassified', pd.NA)
In [62]:
census_4.isnull().sum()
Out[62]:
index                    0
Year                     0
Quarter                  0
NAICS Title           5303
Establishments           0
Wage_million$           28
Employment               0
GDP_trillion$            0
Inflation %              0
Unemployment Rate    25681
LMA                      0
dtype: int64
In [63]:
census_4.apply(lambda x: x.isna().sum()/x.count()*100)
Out[63]:
index                0.000000
Year                 0.000000
Quarter              0.000000
NAICS Title          1.134344
Establishments       0.000000
Wage_million$        0.005923
Employment           0.000000
GDP_trillion$        0.000000
Inflation %          0.000000
Unemployment Rate    5.743687
LMA                  0.000000
dtype: float64
In [64]:
msno.matrix(census_4.sort_values(by='Unemployment Rate'))
Out[64]:
<Axes: >
No description has been provided for this image
In [65]:
census_4 = census_4.dropna(subset=['NAICS Title'])
census_4 = census_4.dropna(subset=['Wage_million$'])
In [66]:
media = census_4['Unemployment Rate'].mean()
census_4['Unemployment Rate'].fillna(media, inplace=True)
In [67]:
census_4.info()

<class 'pandas.core.frame.DataFrame'>
Index: 467467 entries, 0 to 472753
Data columns (total 11 columns):
 #   Column             Non-Null Count   Dtype         
---  ------             --------------   -----         
 0   index              467467 non-null  int64         
 1   Year               467467 non-null  string        
 2   Quarter            467467 non-null  datetime64[ns]
 3   NAICS Title        467467 non-null  category      
 4   Establishments     467467 non-null  int64         
 5   Wage_million$      467467 non-null  float64       
 6   Employment         467467 non-null  int64         
 7   GDP_trillion$      467467 non-null  int64         
 8   Inflation %        467467 non-null  float64       
 9   Unemployment Rate  467467 non-null  float64       
 10  LMA                467467 non-null  category      
dtypes: category(2), datetime64[ns](1), float64(3), int64(4), string(1)
memory usage: 37.0 MB

ANÁLISIS Y GESTIÓN DE OUTLIERS

Para determinar la presencia de outliers en las variables numéricas, en primer lugar, se estudia el coeficiente de asimetría.

Las variables ‘Establishments’, ‘Wage_million’ y ‘Employment’ tienen un coeficiente de asimetría positivo bastante alto, que indica una cola a la derecha de la distribución. También la variable ‘Unemployment’ tiene coeficiente de asimetría positivo pero de forma más contenida. En cambio, la variable ‘GDP_trillion$’ tiene un coeficiente de asimetría cercano a 0 que indica una distribución más o menos simétrica y la variable ‘Inflation %’ tiene un coeficiente de asimetría negativo que indica una cola izquierda en la distribución.

Para gestionar los outliers es necesario separar la variable objetivo del dataset.

Se define la función winzorize with pandas, que tiene como parametro la serie a winsorizar y los percentiles a los cuales ajustar los outliers.

Esta función se incorpora en la función gestiona outliers que incluye las opciones de controlar la incidencia de outliers en el total de datos, winsorizar o convertir los outliers a nulos. En este caso se opta por la opción winsor. Una vez gestionados los outliers, se vuelve a unir el resultado con las variables categóricas y la variable objetivo.

In [68]:
y = census_4['Wage_million$']
census_5 = census_4.drop(['Wage_million$'], axis=1)
In [69]:
census_5.select_dtypes(include=np.number).apply(lambda x: x.skew())
Out[69]:
index                 0.000977
Establishments       12.015200
Employment           15.560795
GDP_trillion$         0.047392
Inflation %          -1.710759
Unemployment Rate     0.918534
dtype: float64
In [70]:
def winsorize_with_pandas(s, limits):
    
    return s.clip(lower=s.quantile(limits[0], interpolation='lower'), 
                  upper=s.quantile(1-limits[1], interpolation='higher'))
In [71]:
def gestiona_outliers(col,clas='winsor'):
    
     print(col.name)
     
     if abs(col.skew()) < 1:
        criterio1 = abs((col-col.mean())/col.std())>3
     else:
        criterio1 = abs((col-col.median()))>8
             
     q1 = col.quantile(0.05)  
       
     q3 = col.quantile(0.95)
     
     IQR=q3-q1
     
     criterio2 = (col<(q1 - 3*IQR))|(col>(q3 + 3*IQR))
     lower = col[criterio1&criterio2&(col<q1)].count()/col.dropna().count()
     upper = col[criterio1&criterio2&(col>q3)].count()/col.dropna().count()
     
     if clas == 'check':
            return(lower*100,upper*100,(lower+upper)*100)
     elif clas == 'winsor':
            return(winsorize_with_pandas(col,(lower,upper)))
     elif clas == 'miss':
            print('\n MissingAntes: ' + str(col.isna().sum()))
            col.loc[criterio1&criterio2] = np.nan
            print('MissingDespues: ' + str(col.isna().sum()) +'\n')
            return(col)

          
census_5_winsor = census_5.select_dtypes(include=np.number).copy().apply(lambda x: gestiona_outliers(x))

index
Establishments
Employment
GDP_trillion$
Inflation %
Unemployment Rate
In [72]:
census_5_winsor.describe()
Out[72]:
index Establishments Employment GDP_trillion$ Inflation % Unemployment Rate
count 467467.000000 467467.000000 467467.000000 4.674670e+05 467467.000000 467467.000000
mean 236360.995936 153.861612 6628.912811 1.296687e+07 1.799484 5.787414
std 136477.703566 428.528981 18385.722831 2.981232e+06 1.426027 1.792872
min 0.000000 1.000000 2.000000 8.412000e+06 -5.600000 2.433333
25% 118150.500000 9.000000 267.000000 1.040547e+07 1.200000 4.433333
50% 236254.000000 26.000000 931.000000 1.260281e+07 2.000000 5.433333
75% 354559.500000 87.000000 3736.000000 1.584316e+07 2.700000 6.900000
max 472762.000000 2892.000000 132014.000000 1.789074e+07 4.500000 14.133333
In [73]:
census_6 = census_5_winsor.join(census_5.select_dtypes(exclude=np.number))
In [74]:
census_6 = census_6.join(y)
In [75]:
var_cont= census_6.select_dtypes(include=np.number)
var_cat= census_6.select_dtypes(exclude=np.number)

INSPECCIÓN GRÁFICA DE LAS VARIABLES CATEGÓRICAS

Observando el histograma de las variables categóricas se nota que en los años 2017 y 2018 hay más observaciones que en los otros años de la serie temporal. Se buscan eventuales valores duplicados y se eliminan.
En cuanto a las variables ‘LMA’ y ‘NAICS Title’ no se observan anomalías en los datos.

In [76]:
for col in var_cat.columns:
    var_cat[col] = var_cat[col].astype(str)
    fig = px.histogram(var_cat, x=col, title=f'Distribución de {col}')
    fig.show()
In [77]:
duplicate_rows = census_4[census_4.duplicated(keep='first')]
duplicate_rows.info()

<class 'pandas.core.frame.DataFrame'>
Index: 0 entries
Data columns (total 11 columns):
 #   Column             Non-Null Count  Dtype         
---  ------             --------------  -----         
 0   index              0 non-null      int64         
 1   Year               0 non-null      string        
 2   Quarter            0 non-null      datetime64[ns]
 3   NAICS Title        0 non-null      category      
 4   Establishments     0 non-null      int64         
 5   Wage_million$      0 non-null      float64       
 6   Employment         0 non-null      int64         
 7   GDP_trillion$      0 non-null      int64         
 8   Inflation %        0 non-null      float64       
 9   Unemployment Rate  0 non-null      float64       
 10  LMA                0 non-null      category      
dtypes: category(2), datetime64[ns](1), float64(3), int64(4), string(1)
memory usage: 11.4 KB
In [78]:
census_6 = census_6.drop_duplicates(keep='first').reset_index(drop=True)

INSPECCIÓN GRÁFICA DE LAS VARIABLES CONTINUAS

Las variables continuas serán agrupadas por ‘Quarter’, dado que cada fecha tiene tantas observaciones como categorías sectoriales.

INFLATION
La gráfica de la tasa de inflación no presenta patrones evidentes. En el intervalo de tiempo observado varía entre en -2 y 4% con un pico de -6% en 2009.

ESTABLISHMENTS
En la gráfica del número de establecimientos se observa un fuerte incremento en 2017. Esto podría ser debido a la mayoría de datos que hay en el conjunto de datos para este período.

WAGE
La gráfica de esta variable deja ver un patrón estacional y tendencia creciente. Normalmente hay un pico de salario en el 4 trimestre. Además, a partir del año 2017 hay un crecimiento anómalo que podría estar correlacionado a la mayor cantidad de datos que hay en el dataset para estos años.

EMPLOYMENT
La gráfica deja ver un patrón estacional en todos los años además de un crecimiento anómalo a partir de 2017 que probablemente esté relacionado con la mayor cantidad de datos que se dispone para este periodo.

GDP
Esta gráfica no presenta patrones estacionales visibles pero presenta una tendencia creciente en los años.

UNEMPLOYMENT RATE
Esta gráfica presenta una tendencia específica pero deja ver un posible patrón estacional: parecería que el 4 trimestre suele ser el con la tasa de desempleo más alta del año.

In [79]:
plt.rcParams["figure.figsize"] = (20, 10)
inflation_byquarter = census_6.groupby('Quarter')['Inflation %'].mean()
inflation_byquarter.plot()
plt.show()
No description has been provided for this image
In [80]:
plt.rcParams["figure.figsize"] = (20, 10)
Establishments_byquarter = census_6.groupby('Quarter')['Establishments'].sum()
Establishments_byquarter.plot()
plt.show()
No description has been provided for this image
In [81]:
plt.rcParams["figure.figsize"] = (20, 10)
Wage_million_byquarter = census_6.groupby('Quarter')['Wage_million$'].sum()
Wage_million_byquarter.plot()
plt.show()
No description has been provided for this image
In [82]:
plt.rcParams["figure.figsize"] = (20, 10)
Employment_byquarter = census_6.groupby('Quarter')['Employment'].sum()
Employment_byquarter.plot()
plt.show()
No description has been provided for this image
In [83]:
plt.rcParams["figure.figsize"] = (20, 10)
GDP_byquarter = census_6.groupby('Quarter')['GDP_trillion$'].mean()
GDP_byquarter.plot()
plt.show()
No description has been provided for this image
In [84]:
plt.rcParams["figure.figsize"] = (20, 10)
Unemployment_byquarter = census_6.groupby('Quarter')['Unemployment Rate'].mean()
Unemployment_byquarter.plot()
plt.show()
No description has been provided for this image

TRASFORMACIÓN DE VARIABLES CATEGÓRICAS A NUMÉRICAS

Con el fin de agregar al modelo también las variables categóricas, se procede a convertirlas en numéricas.

En primer lugar, para la variable categórica ‘LMA’ que tiene 9 categorías se procede a su transformación en dummy a través de la función get dummies.

Para la variable ‘NAICS Title’, que tiene muchas categorías, la gestión en dummies sería complicada, así que se utiliza la función Label Encoder, que permite atribuir un número a cada categoría de la variable.

Una vez que todas las variables sean numéricas se puede estudiar la correlación de cada una con respecto a la variable objetivo.
Las variables con más correlacionadas con la variable objetivo son: ‘Employment’, ‘Establishments’, ‘LMA_New York City Region’.

In [85]:
census_6.head()
Out[85]:
index Establishments Employment GDP_trillion$ Inflation % Unemployment Rate Year Quarter NAICS Title LMA Wage_million$
0 0 6 121 16627803 2.8 6.700000 2018 2018-03-31 Florists Southern Tier Region 0.184196
1 1 30 291 17529646 0.9 5.787414 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting North Country Region 0.667383
2 2 18 591 17529646 0.9 5.787414 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting Capital Region 1.971956
3 3 18 358 17529646 0.9 5.787414 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting Western New York Region 0.930796
4 4 4 65 17529646 0.9 5.787414 2019 2019-03-31 Agriculture, Forestry, Fishing and Hunting New York City Region 0.213331
In [86]:
LMA_encoded = pd.get_dummies(census_6, columns=['LMA'], prefix='LMA', drop_first=True)
census_6 = pd.concat([LMA_encoded],axis=1)
census_6.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 467467 entries, 0 to 467466
Data columns (total 19 columns):
 #   Column                       Non-Null Count   Dtype         
---  ------                       --------------   -----         
 0   index                        467467 non-null  int64         
 1   Establishments               467467 non-null  int64         
 2   Employment                   467467 non-null  int64         
 3   GDP_trillion$                467467 non-null  int64         
 4   Inflation %                  467467 non-null  float64       
 5   Unemployment Rate            467467 non-null  float64       
 6   Year                         467467 non-null  string        
 7   Quarter                      467467 non-null  datetime64[ns]
 8   NAICS Title                  467467 non-null  category      
 9   Wage_million$                467467 non-null  float64       
 10  LMA_Central New York Region  467467 non-null  bool          
 11  LMA_Finger Lakes Region      467467 non-null  bool          
 12  LMA_Hudson Valley Region     467467 non-null  bool          
 13  LMA_Long Island Region       467467 non-null  bool          
 14  LMA_Mohawk Valley Region     467467 non-null  bool          
 15  LMA_New York City Region     467467 non-null  bool          
 16  LMA_North Country Region     467467 non-null  bool          
 17  LMA_Southern Tier Region     467467 non-null  bool          
 18  LMA_Western New York Region  467467 non-null  bool          
dtypes: bool(9), category(1), datetime64[ns](1), float64(3), int64(4), string(1)
memory usage: 37.0 MB
In [87]:
census_6['LMA_Central New York Region'] = census_6['LMA_Central New York Region'].astype(int)
census_6['LMA_Finger Lakes Region'] = census_6['LMA_Finger Lakes Region'].astype(int)
census_6['LMA_Hudson Valley Region'] = census_6['LMA_Hudson Valley Region'].astype(int)
census_6['LMA_Long Island Region'] = census_6['LMA_Long Island Region'].astype(int)
census_6['LMA_Mohawk Valley Region'] = census_6['LMA_Mohawk Valley Region'].astype(int)
census_6['LMA_New York City Region'] = census_6['LMA_New York City Region'].astype(int)
census_6['LMA_North Country Region'] = census_6['LMA_North Country Region'].astype(int)
census_6['LMA_Southern Tier Region'] = census_6['LMA_Southern Tier Region'].astype(int)
census_6['LMA_Western New York Region'] = census_6['LMA_Western New York Region'].astype(int)
In [88]:
le = LabelEncoder()
In [89]:
le.fit(census_6['NAICS Title'].unique())
census_6['NAICS Title']=le.transform(census_6['NAICS Title'])
In [90]:
census_6.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 467467 entries, 0 to 467466
Data columns (total 19 columns):
 #   Column                       Non-Null Count   Dtype         
---  ------                       --------------   -----         
 0   index                        467467 non-null  int64         
 1   Establishments               467467 non-null  int64         
 2   Employment                   467467 non-null  int64         
 3   GDP_trillion$                467467 non-null  int64         
 4   Inflation %                  467467 non-null  float64       
 5   Unemployment Rate            467467 non-null  float64       
 6   Year                         467467 non-null  string        
 7   Quarter                      467467 non-null  datetime64[ns]
 8   NAICS Title                  467467 non-null  int32         
 9   Wage_million$                467467 non-null  float64       
 10  LMA_Central New York Region  467467 non-null  int32         
 11  LMA_Finger Lakes Region      467467 non-null  int32         
 12  LMA_Hudson Valley Region     467467 non-null  int32         
 13  LMA_Long Island Region       467467 non-null  int32         
 14  LMA_Mohawk Valley Region     467467 non-null  int32         
 15  LMA_New York City Region     467467 non-null  int32         
 16  LMA_North Country Region     467467 non-null  int32         
 17  LMA_Southern Tier Region     467467 non-null  int32         
 18  LMA_Western New York Region  467467 non-null  int32         
dtypes: datetime64[ns](1), float64(3), int32(10), int64(4), string(1)
memory usage: 49.9 MB

HACIA LA ESTACIONARIEDAD

Para preparar los datos ante la modelización SARIMAX, es necesario buscar la estacionariedad en la variabele endógena.

Con la estacionariedad significa que las propiedades estadísticas (media, varianza) permanecen constantes a lo largo del tiempo, por lo que si la serie no resulta estacionaria habrá que diferenciarse hasta alcanzar la estacionariedad. Este analisi llevará a la determinación del valor óptimo del parámetro d del modelo SARIMAX.

En primer lugar, se define la función test stationarity, con parámetro la serie temporal, como el test de Dickey-Fuller para evaluar la estacionariedad de la serie temporal.

In [91]:
def test_stationarity(timeseries):
    
    print('Results of Dickey-Fuller Test:')
    dftest = adfuller(timeseries, autolag='AIC', regression ='ct')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    for key,value in dftest[4].items():
        dfoutput['Critical Value (%s)'%key] = value
    print(dfoutput)
In [92]:
test_stationarity(census_6.groupby('Quarter')['Wage_million$'].sum())

Results of Dickey-Fuller Test:
Test Statistic                  0.518118
p-value                         0.996883
#Lags Used                     11.000000
Number of Observations Used    66.000000
Critical Value (1%)            -4.102931
Critical Value (5%)            -3.479243
Critical Value (10%)           -3.167205
dtype: float64

Del test de Dickei-Fuller se observa un p-valor>0.05 que no permite rechazar la H0.

H0: la serie no es estacionaria.

Conclusión: la serie no es estacionaria.

Se procede con la descomposición estacional aplicando la función seasonal decompose para descomponer la serie en componentes de tendencia, estacionales y residuales, utilizando las medias móviles. Esta función permite detectar patrones ocultos de estacionalidad y ayuda a determinar los valores óptimos de los parametros P, D, Q, m del modelo arima.

In [93]:
Wage_million_byquarter_mult = seasonal_decompose(census_4.groupby('Quarter')['Wage_million$'].sum(), model='multiplicative', period=4)
Wage_million_byquarter_mult.plot()
plt.show()
No description has been provided for this image
In [94]:
test_stationarity(Wage_million_byquarter_mult.resid.dropna())

Results of Dickey-Fuller Test:
Test Statistic                -9.832544e+00
p-value                        5.253356e-15
#Lags Used                     1.000000e+00
Number of Observations Used    7.200000e+01
Critical Value (1%)           -4.090351e+00
Critical Value (5%)           -3.473332e+00
Critical Value (10%)          -3.163778e+00
dtype: float64

El test de Dickey-Fuller sobre los residuos devuelve un p-valor < 0, índice de una serie potencialmente estacionaria, como ya se podía apreciar del gráfico de los residuos.

En los próximos pasos se aplicarán unas diferenciaciones a la serie para buscar la estacionariedad.

  1. Estabilización de la varianza a través de una transformación logaritmica.
  2. Diferenciación regular.
  3. Diferenciación estacional.

TRANSFORMACIÓN LOGARITMICA

In [95]:
Wage_million_byquarter_log = np.log(census_6.groupby('Quarter')['Wage_million$'].sum())
Wage_million_byquarter_log.plot()
plt.show()
No description has been provided for this image

Tras la primera diferenciación se observa una menor fluctuación, pero aun no se ha eliminado ni la tendencia ni la componente estacional trimestral que se observaba al principio.
Se procede con la diferenciación regular.

DIFERENCIACIÓN REGULAR

In [96]:
Wage_million_byquarter_logdiff1 = Wage_million_byquarter_log.diff(periods=1)
Wage_million_byquarter_logdiff1.plot()
plt.show()
No description has been provided for this image

En este caso se ha eliminado la tendencia. Sigue siendo evidente el patrón estacional. Para eliminarlo, en el siguiente paso se aplicará una diferenciación estacional.

DIFERENCIACIÓN ESTACIONAL

In [97]:
Wage_million_byquarter_logdiff1_4 = Wage_million_byquarter_logdiff1.diff(periods=4)
Wage_million_byquarter_logdiff1_4.plot()
plt.show()
No description has been provided for this image

Con esta última diferenciación se observa gráficamente que se ha eliminato también la componente estacional.

In [98]:
test_stationarity(Wage_million_byquarter_logdiff1_4.dropna(inplace=False))

Results of Dickey-Fuller Test:
Test Statistic                 -4.664547
p-value                         0.000802
#Lags Used                      8.000000
Number of Observations Used    64.000000
Critical Value (1%)            -4.107677
Critical Value (5%)            -3.481469
Critical Value (10%)           -3.168494
dtype: float64

Tras estas diferenciaciones, el p-valor resultante del test de Dickey-Fuller en la serie diferenciada es < 0.05. Este valor permite rechazar la H0. Ahora la serie es estacionaria y modelizable.

Seguidamente se desarrollará el análisis de autocorrelación que consiste en graficar las funciones de autocorrelación y autocorrelación parcial (ACF y PACF) para identificar posibles relaciones de retardos entre los valores de la serie. Este análisis visual ayuda a determinar los términos autorregresivos (AR) y de media móvil (MA) adecuados (p y q) para el modelo ARIMA, aunque en la fase de modelización se utilizará la herramienta Auto-Arima para obtener el modelo óptimo.

In [99]:
plot_acf(Wage_million_byquarter_logdiff1_4.dropna(inplace=False))
plt.show()

plot_pacf(Wage_million_byquarter_logdiff1_4.dropna(inplace=False), lags=32, method='ywm')
plt.show()
No description has been provided for this image
No description has been provided for this image

Los gráficos de autocorrelación muestran la correlación entre una serie temporal y sus valores pasados. Son una herramienta útil para identificar el orden de un modelo autorregresivo, es decir, los valores pasados que se deben incluir en el modelo.

En lo regular se observan 2 retardos significativos en decrecimiento en ACF y 3 retardos significativos en PACF. En lo estacional se observan 2 retardos significativos en decrecimiento en ACF y 2 retardos significativos en decrecimiento en PACF. Se hipotizan los siguientes ordenes Arima:
(1,1,1)(2,1,2,4)
(1,1,3)(2,1,2,4)

In [100]:
sm.stats.acorr_ljungbox(Wage_million_byquarter_logdiff1_4.dropna(inplace=False), lags=4, return_df=True)
Out[100]:
lb_stat lb_pvalue
1 5.111714 0.023765
2 5.125753 0.077083
3 7.231778 0.064866
4 22.117276 0.000190

En fin, con el test de Ljung-box se obtiene un p-valor > 0.05 que permite no rechazar la H0 en el segundo y tercer retardo.

H0: la serie es incorrelada

En cambio, el primero y cuarto retardo tienen p-valor < 0.05, con lo cual aquí se rechazaría la H0.

PARTICIÓN TRAIN TEST¶

Tras una selección de las variables más significativas, se reparten los datos en el conjunto de entrenamiento y de test: los anteriores a 2017 en el conjunto de train y los sucesivos en el conjunto de test.

El modelo predictivo se basará en los datos trimestrales agrupados, entonces se procede a la agrupación de las variables exógenas y de la endógena sea en el conjunto de train que en el conjunto de test.

In [101]:
census_6['Year'] = pd.to_numeric(census_6['Year'])
In [102]:
train = census_6[census_6['Year'] < 2017].set_index('Quarter').select_dtypes(include=np.number)
train = train.drop('Year', axis=1)
In [103]:
test = census_6[(census_6['Year'] >= 2017)].set_index('Quarter').select_dtypes(include=np.number)
test = test.drop('Year', axis=1)
In [104]:
a=train.groupby('Quarter')[['Establishments', 'Employment']].sum()
b=train.groupby('Quarter')[['Unemployment Rate','Inflation %','GDP_trillion$','LMA_Central New York Region']].mean()
exog_tr=pd.merge(a,b, on='Quarter',how='left')
y_train=train.groupby('Quarter')['Wage_million$'].sum()
tr=pd.concat([exog_tr, y_train], axis=1)
In [105]:
a=test.groupby('Quarter')[['Establishments', 'Employment']].sum()
b=test.groupby('Quarter')[['Unemployment Rate','Inflation %','GDP_trillion$','LMA_Central New York Region']].mean()
exog_tst=pd.merge(a,b, on='Quarter', how='left')
y_test=test.groupby('Quarter')['Wage_million$'].sum()
tst=pd.concat([exog_tst, y_test], axis=1)
In [106]:
exog_tr.shape, exog_tst.shape
Out[106]:
((68, 6), (10, 6))

MODELIZACIÓN

Para representar gráficamente el conjunto de train, test y la predicción se define la función eval model exog que incluye en la representación de la serie temporal, las métricas de evaluación del modelo.

  • El MAPE (Mean Absolute Percentage Error) mide la precisión del modelo, expresando el porcentaje medio de los errores absolutos del modelo con respecto a los valores reales. Un menor valor de la métrica MAPE indíca un modelo más preciso y más semejante al modelo real.

  • El test Ljung-box evalua la presencia de autocorrelación en los residuos, es decir la capacidad de explicabilidad del módelo. Un p-valor > 0.05 permite aceptar la H0.

    H0: los residuos son incorrelados

    Un p-valor bajo indíca la presencia de residuos correlados, que el módelo podría no explicar.

  • R^2 (coeficiente de determinación) es una medida que indica la bondad del ajuste. Cuanto más cercano al 1 esté el coeficiente R^2, mejor el el modelo explica la proporción de varianza en relación con la varianza total de los datos.

En primer lugar se prueban los modelos detectados con el arima manual:
(1,1,1)(2,1,2,4)
(1,1,3)(2,1,2,4)

Sucesivamente se emplea el metodo Auto Arima. Este metodo devuelve la mejor combinación de parametros (p, d, q)(P, D, Q, m), siendo estos parametros:

p: el orden (número de lags temporales) de la parte autorregresiva del modelo.
d: el grado de diferenciación (el número de veces que se han restado los valores consecutivos de la serie).
q: el tamaño de la media móvil del modelo.

P: el orden (número de lags temporales) de la parte estacional del modelo.
D: el grado de diferenciación de la parte estacional del modelo.
Q: el tamaño de la media móvil de la parte estacional del modelo.
m: indica al número de períodos en cada temporada.

El metodo Auto Arima consta de varios parametros, entre los cuales la información acerca de la estacionalidad de la serie y las diferenciaciones regulares (d) y estacionales(D) aplicadas. Contando con estos parametros, el Auto Arima busca automaticamente, entre varias combinaciones posibles de modelos ARIMA, el que minimiza el valor AIC para encontrar el módelo óptimo.

In [107]:
def eval_model_exog(model, tr, tst, exog, name='Model', lags=4):
    lb = np.mean(sm.stats.acorr_ljungbox(model.resid, lags=lags, return_df=True).lb_pvalue)
    
    pred = model.get_forecast(steps=10, exog=exog)
    
    fig1, ax = plt.subplots()
    ax.plot(tr, label='training')
    ax.plot(tst, label='test')
    ax.plot(pred.predicted_mean, label='prediction')
    plt.legend(loc='upper left')
   
    tit = name + ":  LjungBox p-value --> " + str(lb) + "\n MAPE: " + str(round(mean_absolute_percentage_error(tst, pred.predicted_mean)*100, 2)) + "%" + "\n R²: " + str(round(r2_score(tst, pred.predicted_mean), 2))
    
    plt.title(tit)
    plt.ylabel('Wage')
    plt.xlabel('Date')
    plt.show()
    
    print(lb)

MODELO 1 - ARIMA (1,1,1)(2,1,2,4)

In [108]:
model_1= SARIMAX(train.groupby('Quarter')['Wage_million$'].sum(),
                 order=(1,1,1),
                 seasonal_order= (2,1,2,4),
                 exog=exog_tr)

model_1= model_1.fit()
model_1.summary()

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.
Out[108]:
SARIMAX Results
Dep. Variable: Wage_million$ No. Observations: 68
Model: SARIMAX(1, 1, 1)x(2, 1, [1, 2], 4) Log Likelihood -661.112
Date: Tue, 27 Feb 2024 AIC 1348.223
Time: 22:49:42 BIC 1376.084
Sample: 03-31-2000 HQIC 1359.181
- 12-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
Establishments 0.0959 0.310 0.310 0.757 -0.511 0.703
Employment 0.0059 0.007 0.894 0.372 -0.007 0.019
Unemployment Rate -4579.2202 4022.209 -1.138 0.255 -1.25e+04 3304.166
Inflation % -98.6637 959.828 -0.103 0.918 -1979.892 1782.564
GDP_trillion$ 0.0141 0.010 1.386 0.166 -0.006 0.034
LMA_Central New York Region -2.759e+06 3.195 -8.64e+05 0.000 -2.76e+06 -2.76e+06
ar.L1 -0.0300 0.333 -0.090 0.928 -0.683 0.623
ma.L1 -0.5700 0.427 -1.335 0.182 -1.407 0.267
ar.S.L4 0.0686 0.657 0.104 0.917 -1.219 1.356
ar.S.L8 -0.2576 0.473 -0.545 0.586 -1.184 0.669
ma.S.L4 -0.4713 0.650 -0.725 0.469 -1.746 0.803
ma.S.L8 -0.1827 0.582 -0.314 0.754 -1.323 0.958
sigma2 8.988e+07 0.158 5.69e+08 0.000 8.99e+07 8.99e+07
Ljung-Box (L1) (Q): 1.93 Jarque-Bera (JB): 33.87
Prob(Q): 0.16 Prob(JB): 0.00
Heteroskedasticity (H): 0.87 Skew: 1.19
Prob(H) (two-sided): 0.75 Kurtosis: 5.68


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 2.23e+26. Standard errors may be unstable.



In [109]:
eval_model_exog(model_1,y_train,y_test,exog=exog_tst,name='Model_1',lags=4)
No description has been provided for this image 
0.047702985210630225

EVALUACIÓN MODELO 1 En este modelo, el test de Ljung-box no permite aceptar del todo la hipotesis de incorrelación de los residuos, así que el modelo podría no abarcar todos los elementos necesarios para devolver una predicción fiel.
En cuanto a MAPE, el valor cercano al 8% podría ser aceptable.
En fin, el coeficiente de determinación confirma que este modelo es mejor que un no modelo (media de los datos).

MODELO 2 - ARIMA (1,1,3)(2,1,2,4)

In [110]:
model_2= SARIMAX(train.groupby('Quarter')['Wage_million$'].sum(),
                 order=(1,1,3),
                 seasonal_order= (2,1,2,4),
                 exog=exog_tr)

model_2= model_2.fit()
model_2.summary()

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:966: UserWarning:

Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning:

Non-invertible starting MA parameters found. Using zeros as starting parameters.
Out[110]:
SARIMAX Results
Dep. Variable: Wage_million$ No. Observations: 68
Model: SARIMAX(1, 1, 3)x(2, 1, [1, 2], 4) Log Likelihood -655.641
Date: Tue, 27 Feb 2024 AIC 1341.282
Time: 22:49:43 BIC 1373.429
Sample: 03-31-2000 HQIC 1353.926
- 12-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
Establishments 0.0653 0.313 0.209 0.835 -0.548 0.678
Employment 0.0068 0.007 0.982 0.326 -0.007 0.020
Unemployment Rate -4579.2202 4079.232 -1.123 0.262 -1.26e+04 3415.928
Inflation % -98.6636 1073.268 -0.092 0.927 -2202.231 2004.904
GDP_trillion$ 0.0133 0.009 1.418 0.156 -0.005 0.032
LMA_Central New York Region -2.759e+06 3.256 -8.47e+05 0.000 -2.76e+06 -2.76e+06
ar.L1 -0.4920 299.863 -0.002 0.999 -588.213 587.229
ma.L1 -0.3707 299.910 -0.001 0.999 -588.183 587.442
ma.L2 -0.3471 258.740 -0.001 0.999 -507.468 506.774
ma.L3 0.0403 23.305 0.002 0.999 -45.637 45.717
ar.S.L4 0.2921 1.043 0.280 0.779 -1.751 2.336
ar.S.L8 -0.1396 0.721 -0.193 0.847 -1.554 1.274
ma.S.L4 -0.4394 1.079 -0.407 0.684 -2.554 1.675
ma.S.L8 -0.0903 0.886 -0.102 0.919 -1.827 1.647
sigma2 8.852e+07 0.414 2.14e+08 0.000 8.85e+07 8.85e+07
Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 37.71
Prob(Q): 0.88 Prob(JB): 0.00
Heteroskedasticity (H): 1.39 Skew: 0.27
Prob(H) (two-sided): 0.46 Kurtosis: 6.75


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 8e+23. Standard errors may be unstable.



In [111]:
eval_model_exog(model_2,y_train,y_test,exog=exog_tst,name='Model_2',lags=4)
No description has been provided for this image 
0.3078138917708689

Este modelo tiene buena pinta. Se define la función residcheck para evaluar los residuos a través de contrastes de hipotesis.

In [112]:
def residcheck(residuals, lags):

    resid_mean = np.mean(residuals)
    lj_p_val = np.mean(sm.stats.acorr_ljungbox(x=residuals, lags=lags).lb_pvalue)
    norm_p_val =  stats.jarque_bera(residuals)[1]
    adfuller_p = adfuller(residuals)[1]
        
      
    fig = plt.figure(figsize=(10,8))
    layout = (2, 2)
    ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2);
    acf_ax = plt.subplot2grid(layout, (1, 0));
    kde_ax = plt.subplot2grid(layout, (1, 1));
    
    residuals.plot(ax=ts_ax)
    plot_acf(residuals, lags=lags, ax=acf_ax);
    sns.kdeplot(residuals);
    #[ax.set_xlim(1.5) for ax in [acf_ax, kde_ax]]
    sns.despine()
    plt.tight_layout();
    plt.show()
    print("** Mean of the residuals: ", np.around(resid_mean,2))
        
    print("\n** Ljung Box Test, p-value:", np.around(lj_p_val,3), 
        "(>0.05, Uncorrelated)" if (lj_p_val > 0.05) else "(<0.05, Correlated)")
        
    print("\n** Jarque Bera Normality Test, p_value:", np.around(norm_p_val,3),
        "(>0.05, Normal)" if (norm_p_val>0.05) else "(<0.05, Not-normal)")
        
    print("\n** AD Fuller, p_value:", np.around(adfuller_p,3), 
        "(>0.05, Non-stationary)" if (adfuller_p > 0.05) else "(<0.05, Stationary)")
    
    return ts_ax, acf_ax, kde_ax
In [113]:
residcheck(model_2.resid,4)
No description has been provided for this image 
** Mean of the residuals:  -868.83

** Ljung Box Test, p-value: 0.308 (>0.05, Uncorrelated)

** Jarque Bera Normality Test, p_value: 0.0 (<0.05, Not-normal)

** AD Fuller, p_value: 0.0 (<0.05, Stationary)
Out[113]:
(<Axes: xlabel='Quarter'>,
 <Axes: title={'center': 'Autocorrelation'}>,
 <Axes: ylabel='Density'>)

EVALUACIÓN MODELO 2

En este modelo, el test de Ljung-box permite aceptar la hipotesis de incorrelación de los residuos. Así que el modelo tiene una buena explicabilidad y bajo este aspecto permite aceptar con más fuerza la H0 con respecto al modelo 1. En cuanto a MAPE, este modelo tiene un valor ligeramente superior al modelo 1, aunque muy parecido. En fin, también el coeficiente de determinación es muy parecido al modelo 1.

Con respecto al modelo anterior, se prefiere este modelo en cuanto acepta con más fuerza la H0 del test de Ljung-box.

Del analisis de residuos se observa que la media de los residuos es muy lejana a 0. El test Ljung-box confirma incorrelación de los residuos; el test Jarque Bera evidencia una distribución no normal de los residuos; el test Dickey Fuller evidencia estacionariedad de las serie de residuos.

MODELO 3 - AUTO ARIMA (0,1,0)(2,1,2)
Predicción con modelo SARIMAX. En la definición del modelo se necesita indicar el conjunto de train de la variable endogena, el orden arima regular y estacional y el conjunto de entrenamiento de las variables exogenas.

Con el método fit.predict se ajusta el modelo a los datos de entrenamiento y se hace la predicción.

In [114]:
arima_auto = pm.auto_arima(train.groupby('Quarter')[['Wage_million$']].sum(), start_p=1, start_q=1, 
                          test='adf',
                          max_p=5, max_q=5,
                          m=4,
                          d=1,
                          seasonal=True,
                          D=1,
                          trace=True,
                          error_action='ignore',
                          suppress_warnings=True,
                          stepwise=True)

print(arima_auto.summary())

Performing stepwise search to minimize aic
 ARIMA(1,1,1)(1,1,1)[4]             : AIC=1375.049, Time=0.12 sec
 ARIMA(0,1,0)(0,1,0)[4]             : AIC=1374.087, Time=0.01 sec
 ARIMA(1,1,0)(1,1,0)[4]             : AIC=1377.882, Time=0.03 sec
 ARIMA(0,1,1)(0,1,1)[4]             : AIC=1378.575, Time=0.03 sec
 ARIMA(0,1,0)(1,1,0)[4]             : AIC=1375.798, Time=0.02 sec
 ARIMA(0,1,0)(0,1,1)[4]             : AIC=1376.277, Time=0.03 sec
 ARIMA(0,1,0)(1,1,1)[4]             : AIC=1373.164, Time=0.04 sec
 ARIMA(0,1,0)(2,1,1)[4]             : AIC=1369.379, Time=0.07 sec
 ARIMA(0,1,0)(2,1,0)[4]             : AIC=1374.701, Time=0.03 sec
 ARIMA(0,1,0)(2,1,2)[4]             : AIC=1362.910, Time=0.15 sec
 ARIMA(0,1,0)(1,1,2)[4]             : AIC=1369.655, Time=0.09 sec
 ARIMA(1,1,0)(2,1,2)[4]             : AIC=1365.342, Time=0.20 sec
 ARIMA(0,1,1)(2,1,2)[4]             : AIC=1364.928, Time=0.12 sec
 ARIMA(1,1,1)(2,1,2)[4]             : AIC=1366.908, Time=0.45 sec
 ARIMA(0,1,0)(2,1,2)[4] intercept   : AIC=1368.003, Time=0.11 sec

Best model:  ARIMA(0,1,0)(2,1,2)[4]          
Total fit time: 1.505 seconds
                                       SARIMAX Results                                        
==============================================================================================
Dep. Variable:                                      y   No. Observations:                   68
Model:             SARIMAX(0, 1, 0)x(2, 1, [1, 2], 4)   Log Likelihood                -676.455
Date:                                Tue, 27 Feb 2024   AIC                           1362.910
Time:                                        22:49:45   BIC                           1373.626
Sample:                                    03-31-2000   HQIC                          1367.125
                                         - 12-31-2016                                         
Covariance Type:                                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.S.L4        1.1932      0.129      9.217      0.000       0.939       1.447
ar.S.L8       -0.7101      0.075     -9.420      0.000      -0.858      -0.562
ma.S.L4       -1.5177      0.221     -6.875      0.000      -1.950      -1.085
ma.S.L8        0.8349      0.194      4.303      0.000       0.455       1.215
sigma2      1.228e+08   2.07e-09   5.93e+16      0.000    1.23e+08    1.23e+08
===================================================================================
Ljung-Box (L1) (Q):                   7.19   Jarque-Bera (JB):                 2.26
Prob(Q):                              0.01   Prob(JB):                         0.32
Heteroskedasticity (H):               1.30   Skew:                            -0.35
Prob(H) (two-sided):                  0.56   Kurtosis:                         3.62
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 2.48e+33. Standard errors may be unstable.

El modelo óptimo devuelto por Auto Arima es (0,1,0) (2,1,2,4).

In [115]:
model_3= SARIMAX(train.groupby('Quarter')['Wage_million$'].sum(),
                 order=(0,1,0),
                 seasonal_order= (2,1,2,4),
                 exog=exog_tr)

model_3= model_3.fit()
model_3.summary()

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning:

Maximum Likelihood optimization failed to converge. Check mle_retvals
Out[115]:
SARIMAX Results
Dep. Variable: Wage_million$ No. Observations: 68
Model: SARIMAX(0, 1, 0)x(2, 1, [1, 2], 4) Log Likelihood -669.530
Date: Tue, 27 Feb 2024 AIC 1361.060
Time: 22:49:46 BIC 1384.635
Sample: 03-31-2000 HQIC 1370.332
- 12-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
Establishments -0.1298 0.251 -0.517 0.605 -0.622 0.363
Employment 0.0098 0.007 1.325 0.185 -0.005 0.024
Unemployment Rate -4579.2202 4163.376 -1.100 0.271 -1.27e+04 3580.848
Inflation % -98.6637 703.405 -0.140 0.888 -1477.312 1279.985
GDP_trillion$ 0.0203 0.009 2.362 0.018 0.003 0.037
LMA_Central New York Region -2.759e+06 3.775 -7.31e+05 0.000 -2.76e+06 -2.76e+06
ar.S.L4 0.3414 0.493 0.693 0.488 -0.624 1.307
ar.S.L8 -0.3310 0.462 -0.717 0.474 -1.236 0.574
ma.S.L4 -0.3943 0.474 -0.831 0.406 -1.324 0.536
ma.S.L8 -0.0204 0.451 -0.045 0.964 -0.904 0.863
sigma2 1.044e+08 0.284 3.68e+08 0.000 1.04e+08 1.04e+08
Ljung-Box (L1) (Q): 13.76 Jarque-Bera (JB): 2.35
Prob(Q): 0.00 Prob(JB): 0.31
Heteroskedasticity (H): 1.90 Skew: -0.23
Prob(H) (two-sided): 0.15 Kurtosis: 3.82


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 3.31e+24. Standard errors may be unstable.



In [116]:
eval_model_exog(model_3,y_train,y_test,exog=exog_tst,name='Model_3',lags=4)
No description has been provided for this image 
0.0051786510385925365

EVALUACIÓN MODELO 3

Este modelo presenta un test de Ljung-box que no permite aceptar la H0 de incorrelación de la serie. En cambio, el MAPE es aceptable, ya que indica una media de errores absolutos del 6,03%. El coeficiente de determinación indica que el modelo da predicciones mejores al no modelo o media de los datos.

En los próximos pasos se intentará mejorar este modelo aportado una diferenciación adicional en la parte regular y estacional.

MODELO 4 - (2,2,2)(2,1,2)
Se aporta una diferenciación adicional en lo regular para mejorar el test Ljung-Box con respecto al modelo anterior (modelo 3).

In [117]:
arima_auto = pm.auto_arima(train.groupby('Quarter')[['Wage_million$']].sum(), start_p=1, start_q=1,
                          test='adf',
                          max_p=5, max_q=5,
                          m=4,
                          d=2,
                          seasonal=True,
                          D=1,
                          trace=True,
                          error_action='ignore',
                          suppress_warnings=True,
                          stepwise=True)

print(arima_auto.summary())

Performing stepwise search to minimize aic
 ARIMA(1,2,1)(1,1,1)[4]             : AIC=1351.132, Time=0.08 sec
 ARIMA(0,2,0)(0,1,0)[4]             : AIC=1417.088, Time=0.01 sec
 ARIMA(1,2,0)(1,1,0)[4]             : AIC=1409.451, Time=0.03 sec
 ARIMA(0,2,1)(0,1,1)[4]             : AIC=inf, Time=0.05 sec
 ARIMA(1,2,1)(0,1,1)[4]             : AIC=1351.911, Time=0.09 sec
 ARIMA(1,2,1)(1,1,0)[4]             : AIC=1351.925, Time=0.05 sec
 ARIMA(1,2,1)(2,1,1)[4]             : AIC=1350.510, Time=0.14 sec
 ARIMA(1,2,1)(2,1,0)[4]             : AIC=1351.790, Time=0.08 sec
 ARIMA(1,2,1)(2,1,2)[4]             : AIC=1348.978, Time=0.17 sec
 ARIMA(1,2,1)(1,1,2)[4]             : AIC=1350.905, Time=0.11 sec
 ARIMA(0,2,1)(2,1,2)[4]             : AIC=1355.226, Time=0.18 sec
 ARIMA(1,2,0)(2,1,2)[4]             : AIC=1394.110, Time=0.22 sec
 ARIMA(2,2,1)(2,1,2)[4]             : AIC=1349.966, Time=0.35 sec
 ARIMA(1,2,2)(2,1,2)[4]             : AIC=1349.775, Time=0.28 sec
 ARIMA(0,2,0)(2,1,2)[4]             : AIC=1403.573, Time=0.10 sec
 ARIMA(0,2,2)(2,1,2)[4]             : AIC=1351.135, Time=0.29 sec
 ARIMA(2,2,0)(2,1,2)[4]             : AIC=1380.568, Time=0.17 sec
 ARIMA(2,2,2)(2,1,2)[4]             : AIC=1345.791, Time=0.44 sec
 ARIMA(2,2,2)(1,1,2)[4]             : AIC=1348.432, Time=0.29 sec
 ARIMA(2,2,2)(2,1,1)[4]             : AIC=1347.977, Time=0.28 sec
 ARIMA(2,2,2)(1,1,1)[4]             : AIC=1349.104, Time=0.16 sec
 ARIMA(3,2,2)(2,1,2)[4]             : AIC=inf, Time=0.48 sec
 ARIMA(2,2,3)(2,1,2)[4]             : AIC=inf, Time=0.49 sec
 ARIMA(1,2,3)(2,1,2)[4]             : AIC=1354.613, Time=0.39 sec
 ARIMA(3,2,1)(2,1,2)[4]             : AIC=1351.966, Time=0.32 sec
 ARIMA(3,2,3)(2,1,2)[4]             : AIC=inf, Time=0.58 sec
 ARIMA(2,2,2)(2,1,2)[4] intercept   : AIC=1348.945, Time=0.59 sec

Best model:  ARIMA(2,2,2)(2,1,2)[4]          
Total fit time: 6.446 seconds
                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                                 y   No. Observations:                   68
Model:             SARIMAX(2, 2, 2)x(2, 1, 2, 4)   Log Likelihood                -663.896
Date:                           Tue, 27 Feb 2024   AIC                           1345.791
Time:                                   22:49:53   BIC                           1364.935
Sample:                               03-31-2000   HQIC                          1353.308
                                    - 12-31-2016                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -1.3030      0.249     -5.239      0.000      -1.790      -0.815
ar.L2         -0.4033      0.194     -2.081      0.037      -0.783      -0.024
ma.L1          0.0066      0.493      0.013      0.989      -0.960       0.973
ma.L2         -0.9594      0.400     -2.396      0.017      -1.744      -0.175
ar.S.L4        1.2160      0.282      4.313      0.000       0.663       1.769
ar.S.L8       -0.7520      0.228     -3.304      0.001      -1.198      -0.306
ma.S.L4       -1.4020      0.376     -3.733      0.000      -2.138      -0.666
ma.S.L8        0.7950      0.341      2.329      0.020       0.126       1.464
sigma2       1.63e+08   3.87e-09   4.21e+16      0.000    1.63e+08    1.63e+08
===================================================================================
Ljung-Box (L1) (Q):                   0.56   Jarque-Bera (JB):                 2.16
Prob(Q):                              0.45   Prob(JB):                         0.34
Heteroskedasticity (H):               1.07   Skew:                            -0.08
Prob(H) (two-sided):                  0.88   Kurtosis:                         3.90
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 1.32e+33. Standard errors may be unstable.
In [118]:
model_4= SARIMAX(train.groupby('Quarter')['Wage_million$'].sum(),
                 order=(2,2,2),
                 seasonal_order= (2,1,2,4),
                 exog=exog_tr)

model_4= model_4.fit()
model_4.summary()

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:966: UserWarning:

Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning:

Non-invertible starting MA parameters found. Using zeros as starting parameters.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning:

Maximum Likelihood optimization failed to converge. Check mle_retvals
Out[118]:
SARIMAX Results
Dep. Variable: Wage_million$ No. Observations: 68
Model: SARIMAX(2, 2, 2)x(2, 1, 2, 4) Log Likelihood -663.612
Date: Tue, 27 Feb 2024 AIC 1357.224
Time: 22:49:53 BIC 1389.131
Sample: 03-31-2000 HQIC 1369.752
- 12-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
Establishments -0.0277 0.514 -0.054 0.957 -1.035 0.980
Employment 0.0127 0.012 1.068 0.286 -0.011 0.036
Unemployment Rate -448.5351 6512.963 -0.069 0.945 -1.32e+04 1.23e+04
Inflation % -795.6796 998.332 -0.797 0.425 -2752.374 1161.015
GDP_trillion$ 0.0216 0.012 1.754 0.079 -0.003 0.046
LMA_Central New York Region -5.541e+06 9.112 -6.08e+05 0.000 -5.54e+06 -5.54e+06
ar.L1 -1.2160 0.448 -2.714 0.007 -2.094 -0.338
ar.L2 -0.4173 0.307 -1.361 0.173 -1.018 0.184
ma.L1 -0.1237 0.667 -0.186 0.853 -1.430 1.183
ma.L2 -0.8639 0.512 -1.688 0.091 -1.867 0.139
ar.S.L4 0.4203 9.246 0.045 0.964 -17.702 18.543
ar.S.L8 0.4596 4.589 0.100 0.920 -8.534 9.453
ma.S.L4 -0.3076 8.918 -0.034 0.972 -17.786 17.171
ma.S.L8 -0.5724 5.322 -0.108 0.914 -11.004 9.859
sigma2 1.575e+08 0.023 6.91e+09 0.000 1.58e+08 1.58e+08
Ljung-Box (L1) (Q): 0.01 Jarque-Bera (JB): 10.21
Prob(Q): 0.91 Prob(JB): 0.01
Heteroskedasticity (H): 2.46 Skew: -0.14
Prob(H) (two-sided): 0.04 Kurtosis: 4.97


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 8.6e+26. Standard errors may be unstable.



In [119]:
eval_model_exog(model_4,y_train,y_test,exog=exog_tst,name='Model_4',lags=4)
No description has been provided for this image 
0.2270505892206734

EVALUACIÓN MODELO 4

Añadiendo una diferenciación regular se ha mejorado la métrica Ljung-box, permitiendo ahora aceptar la H0 de incorrelación de los resiudos. En cambio se ha empeorado el MAPE, que ha subido a casi el 23%. Además, el coeficiente de determinación indica que este modelo es peor que un no modelo osea la media de los datos.

MODELO 5 - SARIMAX CON UNA DIFERENCIACIÓN REGULAR Y DOS ESTACIONALES
En este modelo se añade una diferenciación estacional más con respecto al modelo 3.

In [120]:
arima_auto = pm.auto_arima(train.groupby('Quarter')[['Wage_million$']].sum(), start_p=1, start_q=1,
                          test='adf',
                          max_p=5, max_q=5,
                          m=4,
                          d=1,
                          seasonal=True,
                          D=2,
                          trace=True,
                          error_action='ignore',
                          suppress_warnings=True,
                          stepwise=True)

print(arima_auto.summary())

Performing stepwise search to minimize aic
 ARIMA(1,1,1)(1,2,1)[4]             : AIC=1301.521, Time=0.08 sec
 ARIMA(0,1,0)(0,2,0)[4]             : AIC=1335.732, Time=0.01 sec
 ARIMA(1,1,0)(1,2,0)[4]             : AIC=1333.908, Time=0.03 sec
 ARIMA(0,1,1)(0,2,1)[4]             : AIC=inf, Time=0.09 sec
 ARIMA(1,1,1)(0,2,1)[4]             : AIC=1299.543, Time=0.04 sec
 ARIMA(1,1,1)(0,2,0)[4]             : AIC=1338.266, Time=0.03 sec
 ARIMA(1,1,1)(0,2,2)[4]             : AIC=1301.471, Time=0.08 sec
 ARIMA(1,1,1)(1,2,0)[4]             : AIC=1335.786, Time=0.09 sec
 ARIMA(1,1,1)(1,2,2)[4]             : AIC=inf, Time=0.23 sec
 ARIMA(1,1,0)(0,2,1)[4]             : AIC=inf, Time=0.10 sec
 ARIMA(2,1,1)(0,2,1)[4]             : AIC=1300.361, Time=0.18 sec
 ARIMA(1,1,2)(0,2,1)[4]             : AIC=inf, Time=0.24 sec
 ARIMA(0,1,0)(0,2,1)[4]             : AIC=inf, Time=0.06 sec
 ARIMA(0,1,2)(0,2,1)[4]             : AIC=inf, Time=0.08 sec
 ARIMA(2,1,0)(0,2,1)[4]             : AIC=1299.569, Time=0.05 sec
 ARIMA(2,1,2)(0,2,1)[4]             : AIC=1303.854, Time=0.14 sec
 ARIMA(1,1,1)(0,2,1)[4] intercept   : AIC=1301.164, Time=0.14 sec

Best model:  ARIMA(1,1,1)(0,2,1)[4]          
Total fit time: 1.708 seconds
                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                                 y   No. Observations:                   68
Model:             SARIMAX(1, 1, 1)x(0, 2, 1, 4)   Log Likelihood                -645.771
Date:                           Tue, 27 Feb 2024   AIC                           1299.543
Time:                                   22:49:55   BIC                           1307.853
Sample:                               03-31-2000   HQIC                          1302.787
                                    - 12-31-2016                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.0074      0.629      0.012      0.991      -1.225       1.239
ma.L1         -0.5836      0.506     -1.152      0.249      -1.576       0.409
ma.S.L4       -0.8956      0.237     -3.784      0.000      -1.360      -0.432
sigma2      3.054e+08   1.27e-09   2.41e+17      0.000    3.05e+08    3.05e+08
===================================================================================
Ljung-Box (L1) (Q):                   0.00   Jarque-Bera (JB):               185.04
Prob(Q):                              0.98   Prob(JB):                         0.00
Heteroskedasticity (H):               1.27   Skew:                            -1.72
Prob(H) (two-sided):                  0.59   Kurtosis:                        10.96
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 2.64e+33. Standard errors may be unstable.
In [121]:
model_5= SARIMAX(train.groupby('Quarter')['Wage_million$'].sum(),
                 order=(0,1,1),
                 seasonal_order= (0,2,1,4),
                 exog=exog_tr)

model_5= model_5.fit()
model_5.summary()

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.
Out[121]:
SARIMAX Results
Dep. Variable: Wage_million$ No. Observations: 68
Model: SARIMAX(0, 1, 1)x(0, 2, 1, 4) Log Likelihood -631.951
Date: Tue, 27 Feb 2024 AIC 1281.902
Time: 22:49:56 BIC 1300.600
Sample: 03-31-2000 HQIC 1289.201
- 12-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
Establishments 0.0501 0.524 0.096 0.924 -0.977 1.078
Employment 0.0052 0.012 0.426 0.670 -0.019 0.029
Unemployment Rate -5507.5457 6029.413 -0.913 0.361 -1.73e+04 6309.886
Inflation % -123.6356 1599.727 -0.077 0.938 -3259.042 3011.771
GDP_trillion$ 0.0135 0.017 0.812 0.417 -0.019 0.046
LMA_Central New York Region -2.135e+06 1.162 -1.84e+06 0.000 -2.14e+06 -2.14e+06
ma.L1 -0.8786 0.443 -1.983 0.047 -1.747 -0.010
ma.S.L4 -0.8432 0.503 -1.676 0.094 -1.829 0.143
sigma2 2.021e+08 0.183 1.11e+09 0.000 2.02e+08 2.02e+08
Ljung-Box (L1) (Q): 0.12 Jarque-Bera (JB): 113.65
Prob(Q): 0.73 Prob(JB): 0.00
Heteroskedasticity (H): 0.92 Skew: -1.56
Prob(H) (two-sided): 0.85 Kurtosis: 9.04


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 9.89e+24. Standard errors may be unstable.



In [122]:
eval_model_exog(model_5,y_train,y_test,exog=exog_tst,name='Model_5',lags=4)
No description has been provided for this image 
0.7209935005413565
In [123]:
residcheck(model_5.resid,4)
No description has been provided for this image 
** Mean of the residuals:  -1468.15

** Ljung Box Test, p-value: 0.721 (>0.05, Uncorrelated)

** Jarque Bera Normality Test, p_value: 0.0 (<0.05, Not-normal)

** AD Fuller, p_value: 0.0 (<0.05, Stationary)
Out[123]:
(<Axes: xlabel='Quarter'>,
 <Axes: title={'center': 'Autocorrelation'}>,
 <Axes: ylabel='Density'>)

EVALUACIÓN MODELO 5

En el modelo 5, el test de Ljung-box permite aceptar con fuerza que los residuos son incorrelados. El error medio cercano al 6% es aceptable y el coeficiente de determinación indica que el modelo tiene un buen ajuste.

Del analisis de residuos se observa que la media de los residuos es muy lejana a 0. El test Ljung-box confirma incorrelación de los residuos; el test Jarque Bera evidencia una distribución no normal de los residuos; el test Dickey Fuller evidencia estacionariedad de las serie de residuos.

MODELO 6 - SARIMAX CON DOBLE DIFERENCIACIÓN EN LO REGULAR Y ESTACIONAL

En este modelo se añaden una deferenciación más sea en lo regular que en lo estacional, con respecto al modelo 3.

In [124]:
arima_auto = pm.auto_arima(train.groupby('Quarter')[['Wage_million$']].sum(), start_p=1, start_q=1,
                          test='adf',
                          max_p=5, max_q=5,
                          m=4,
                          d=2,
                          seasonal=True,
                          D=2,
                          trace=True,
                          error_action='ignore',
                          suppress_warnings=True,
                          stepwise=True)

print(arima_auto.summary())

Performing stepwise search to minimize aic
 ARIMA(1,2,1)(1,2,1)[4]             : AIC=1305.738, Time=0.07 sec
 ARIMA(0,2,0)(0,2,0)[4]             : AIC=1374.286, Time=0.01 sec
 ARIMA(1,2,0)(1,2,0)[4]             : AIC=1345.844, Time=0.03 sec
 ARIMA(0,2,1)(0,2,1)[4]             : AIC=inf, Time=0.10 sec
 ARIMA(1,2,1)(0,2,1)[4]             : AIC=1304.770, Time=0.04 sec
 ARIMA(1,2,1)(0,2,0)[4]             : AIC=1311.904, Time=0.04 sec
 ARIMA(1,2,1)(0,2,2)[4]             : AIC=1304.920, Time=0.13 sec
 ARIMA(1,2,1)(1,2,0)[4]             : AIC=1313.305, Time=0.08 sec
 ARIMA(1,2,1)(1,2,2)[4]             : AIC=1306.912, Time=0.18 sec
 ARIMA(1,2,0)(0,2,1)[4]             : AIC=inf, Time=0.14 sec
 ARIMA(2,2,1)(0,2,1)[4]             : AIC=1299.347, Time=0.06 sec
 ARIMA(2,2,1)(0,2,0)[4]             : AIC=1310.152, Time=0.13 sec
 ARIMA(2,2,1)(1,2,1)[4]             : AIC=1300.991, Time=0.12 sec
 ARIMA(2,2,1)(0,2,2)[4]             : AIC=1300.500, Time=0.15 sec
 ARIMA(2,2,1)(1,2,0)[4]             : AIC=1311.384, Time=0.13 sec
 ARIMA(2,2,1)(1,2,2)[4]             : AIC=1302.470, Time=0.26 sec
 ARIMA(2,2,0)(0,2,1)[4]             : AIC=1301.755, Time=0.04 sec
 ARIMA(3,2,1)(0,2,1)[4]             : AIC=1300.885, Time=0.12 sec
 ARIMA(2,2,2)(0,2,1)[4]             : AIC=1296.946, Time=0.13 sec
 ARIMA(2,2,2)(0,2,0)[4]             : AIC=1300.594, Time=0.08 sec
 ARIMA(2,2,2)(1,2,1)[4]             : AIC=1296.735, Time=0.14 sec
 ARIMA(2,2,2)(1,2,0)[4]             : AIC=1302.241, Time=0.18 sec
 ARIMA(2,2,2)(2,2,1)[4]             : AIC=1297.169, Time=0.23 sec
 ARIMA(2,2,2)(1,2,2)[4]             : AIC=1297.838, Time=0.30 sec
 ARIMA(2,2,2)(0,2,2)[4]             : AIC=1295.916, Time=0.15 sec
 ARIMA(1,2,2)(0,2,2)[4]             : AIC=1305.723, Time=0.20 sec
 ARIMA(3,2,2)(0,2,2)[4]             : AIC=1296.673, Time=0.27 sec
 ARIMA(2,2,3)(0,2,2)[4]             : AIC=1300.367, Time=0.23 sec
 ARIMA(1,2,3)(0,2,2)[4]             : AIC=1311.200, Time=0.24 sec
 ARIMA(3,2,1)(0,2,2)[4]             : AIC=1301.449, Time=0.22 sec
 ARIMA(3,2,3)(0,2,2)[4]             : AIC=1302.104, Time=0.45 sec
 ARIMA(2,2,2)(0,2,2)[4] intercept   : AIC=1299.056, Time=0.15 sec

Best model:  ARIMA(2,2,2)(0,2,2)[4]          
Total fit time: 4.812 seconds
                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                                 y   No. Observations:                   68
Model:             SARIMAX(2, 2, 2)x(0, 2, 2, 4)   Log Likelihood                -640.958
Date:                           Tue, 27 Feb 2024   AIC                           1295.916
Time:                                   22:50:01   BIC                           1310.339
Sample:                               03-31-2000   HQIC                          1301.534
                                    - 12-31-2016                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -1.2998      0.306     -4.251      0.000      -1.899      -0.701
ar.L2         -0.5029      0.259     -1.945      0.052      -1.010       0.004
ma.L1          0.0805      0.416      0.194      0.846      -0.734       0.895
ma.L2         -0.8016      0.313     -2.561      0.010      -1.415      -0.188
ma.S.L4       -0.6090      0.393     -1.550      0.121      -1.379       0.161
ma.S.L8       -0.2189      0.208     -1.051      0.293      -0.627       0.189
sigma2       3.73e+08   6.99e-10   5.33e+17      0.000    3.73e+08    3.73e+08
===================================================================================
Ljung-Box (L1) (Q):                   0.29   Jarque-Bera (JB):                 7.39
Prob(Q):                              0.59   Prob(JB):                         0.02
Heteroskedasticity (H):               0.85   Skew:                            -0.35
Prob(H) (two-sided):                  0.73   Kurtosis:                         4.60
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 7.02e+33. Standard errors may be unstable.
In [125]:
model_6= SARIMAX(train.groupby('Quarter')['Wage_million$'].sum(),
                 order=(2,2,2),
                 seasonal_order= (0,2,2,4),
                 exog=exog_tr)

model_6= model_6.fit()
model_6.summary()

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:1009: UserWarning:

Non-invertible starting seasonal moving average Using zeros as starting parameters.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning:

Maximum Likelihood optimization failed to converge. Check mle_retvals
Out[125]:
SARIMAX Results
Dep. Variable: Wage_million$ No. Observations: 68
Model: SARIMAX(2, 2, 2)x(0, 2, 2, 4) Log Likelihood -640.631
Date: Tue, 27 Feb 2024 AIC 1307.262
Time: 22:50:02 BIC 1334.048
Sample: 03-31-2000 HQIC 1317.696
- 12-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
Establishments 0.0757 1.116 0.068 0.946 -2.111 2.262
Employment 0.0086 0.026 0.339 0.735 -0.041 0.059
Unemployment Rate -1829.9405 1.18e+04 -0.155 0.877 -2.49e+04 2.13e+04
Inflation % -869.8040 2068.617 -0.420 0.674 -4924.218 3184.610
GDP_trillion$ 0.0248 0.027 0.921 0.357 -0.028 0.078
LMA_Central New York Region -5.268e+06 21.994 -2.4e+05 0.000 -5.27e+06 -5.27e+06
ar.L1 -1.0465 3.008 -0.348 0.728 -6.941 4.848
ar.L2 -0.3235 0.942 -0.343 0.731 -2.170 1.523
ma.L1 -0.2302 3.149 -0.073 0.942 -6.402 5.942
ma.L2 -0.6610 3.030 -0.218 0.827 -6.599 5.277
ma.S.L4 -0.7734 0.931 -0.831 0.406 -2.598 1.051
ma.S.L8 -0.0328 0.435 -0.075 0.940 -0.885 0.819
sigma2 4.264e+08 0.062 6.87e+09 0.000 4.26e+08 4.26e+08
Ljung-Box (L1) (Q): 0.37 Jarque-Bera (JB): 21.24
Prob(Q): 0.54 Prob(JB): 0.00
Heteroskedasticity (H): 1.45 Skew: -0.76
Prob(H) (two-sided): 0.42 Kurtosis: 5.54


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 7.24e+26. Standard errors may be unstable.



In [126]:
eval_model_exog(model_6,y_train,y_test,exog=exog_tst,name='Model_6',lags=4)
No description has been provided for this image 
0.193171780731735

EVALUACIÓN MODELO 6
En este modelo, el test de Ljung-box permite aceptar la hipotesis de incorrelación de los residuos. En cuanto a MAPE, este modelo tiene un valor alto aprox. 14% y el coeficiente de determinación indica que este modelo no es mejor de un no modelo o media de los datos.

MODELO 7 - (0,1,0)(2,1,2,4) SIN VARIABLES EXÓGENAS

Con este modelo se hace una predicción de la serie temporal sin considerar las variables exógenas y considerando los parametros óptimos del modelo Arima con una sola diferenciación sea en la parte regular que en la parte estacional. (0,1,0)(2,1,2,4)

Se define también la función eval model end para la visualización gráfica del modelo predictivo.

In [127]:
def eval_model_end(model,tr,tst,name='Model',lags=4):
    lb = np.mean(sm.stats.acorr_ljungbox(model.resid, lags=lags, return_df=True).lb_pvalue)
    pred = model.forecast(steps=len(tst))
    fig1, ax = plt.subplots()
    ax.plot(tr, label='training')
    ax.plot(tst, label='test')
    ax.plot(pred, label='prediction')
    plt.legend(loc='upper left')
    tit = (name + ":  LjungBox p-value --> " + str(lb) + "\n MAPE: " + str(round(mean_absolute_percentage_error(tst, pred) * 100, 2)) + "%" + "\n R²: " + str(round(r2_score(tst, pred), 2)))    
    plt.title(tit)
    plt.ylabel('Pasajeros')
    plt.xlabel('Date')
    plt.show()
    print(lb)
In [128]:
model_7= SARIMAX(train.groupby('Quarter')['Wage_million$'].sum(),
                 order=(0,1,0),
                 seasonal_order= (2,1,2,4))

model_7= model_7.fit()
model_7.summary()

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.

C:\Users\giuli\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

No frequency information was provided, so inferred frequency Q-DEC will be used.
Out[128]:
SARIMAX Results
Dep. Variable: Wage_million$ No. Observations: 68
Model: SARIMAX(0, 1, 0)x(2, 1, [1, 2], 4) Log Likelihood -676.455
Date: Tue, 27 Feb 2024 AIC 1362.910
Time: 22:50:03 BIC 1373.626
Sample: 03-31-2000 HQIC 1367.125
- 12-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.S.L4 1.1932 0.129 9.217 0.000 0.939 1.447
ar.S.L8 -0.7101 0.075 -9.420 0.000 -0.858 -0.562
ma.S.L4 -1.5177 0.221 -6.875 0.000 -1.950 -1.085
ma.S.L8 0.8349 0.194 4.303 0.000 0.455 1.215
sigma2 1.228e+08 2.07e-09 5.93e+16 0.000 1.23e+08 1.23e+08
Ljung-Box (L1) (Q): 7.19 Jarque-Bera (JB): 2.26
Prob(Q): 0.01 Prob(JB): 0.32
Heteroskedasticity (H): 1.30 Skew: -0.35
Prob(H) (two-sided): 0.56 Kurtosis: 3.62


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 2.48e+33. Standard errors may be unstable.



In [129]:
eval_model_end(model_7,y_train,y_test,name='Model_7',lags=4)
No description has been provided for this image 
0.23908083913032324

EVALUACIÓN MODELO 7

El modelo 7 es bastante peor de los anteriores modelos que empleaban también las variables exógenas. Aunque se puede apreciar del test de Ljung-box que el módelo tiene una buena explicabilidad de los residuos, ya que se acepta la H0 la cual indica la incorrelación de los mismos, el MAPE es muy elevado. Este valor de MAPE es sintoma de una predicción poco precisa, como se puede apreciar también de la gráfica. La presencia de las variables exógenas es muy importante para afinar la predicción, ya que sin ellas se habría obtenido un modelo no utilizable.

MODELO 8 - RANDOM FOREST

En primer lugar se adapta el conjunto de datos incluyendo las columnas 'Year' y 'Month' para aportar una referencia temporal al modelo.

Se crea un objeto de modelo de regresión utilizando el algoritmo Random Forest y se adapta a los datos de entrenamiento.
Seguidamente, con la libreria GridSearchCV se determinan los parametros óptimos para la creación del modelo según la métrica neg_mean_squared_error. Estos parametros se incluyen en el nuevo modelo best model y tras haber hecho la predicción del conjunto de test, se representan gráficamente.

In [130]:
df = pd.concat([tr, tst], axis=0)
df.head()
Out[130]:
Establishments Employment Unemployment Rate Inflation % GDP_trillion$ LMA_Central New York Region Wage_million$
Quarter
2000-03-31 828153 35827751 5.215549 3.4 8412000.0 0.080363 187630.370705
2000-06-30 828346 36744618 3.997643 1.8 8412000.0 0.080798 157076.699272
2000-09-30 831094 36704144 3.879867 2.5 8412000.0 0.080814 155927.064461
2000-12-31 825910 37161011 3.844675 2.2 8412000.0 0.080782 183365.782351
2001-03-31 829990 36111104 4.804509 2.7 8783000.0 0.081250 204286.990668
In [131]:
df['Year']=df.index.astype(str).str[:4].astype(int)
df['Month']=df.index.astype(str).str[5:7].astype(int)
In [132]:
df.head()
Out[132]:
Establishments Employment Unemployment Rate Inflation % GDP_trillion$ LMA_Central New York Region Wage_million$ Year Month
Quarter
2000-03-31 828153 35827751 5.215549 3.4 8412000.0 0.080363 187630.370705 2000 3
2000-06-30 828346 36744618 3.997643 1.8 8412000.0 0.080798 157076.699272 2000 6
2000-09-30 831094 36704144 3.879867 2.5 8412000.0 0.080814 155927.064461 2000 9
2000-12-31 825910 37161011 3.844675 2.2 8412000.0 0.080782 183365.782351 2000 12
2001-03-31 829990 36111104 4.804509 2.7 8783000.0 0.081250 204286.990668 2001 3
In [133]:
df = df.sort_values(by='Quarter')
train_ml = df.iloc[:68]
test_ml = df.iloc[-10:]
y_train_ml = train_ml['Wage_million$']
y_test_ml = test_ml['Wage_million$']
exog_tr_ml = train_ml.drop(columns=['Wage_million$'])
exog_tst_ml = test_ml.drop(columns=['Wage_million$'])
In [134]:
model_8 = RandomForestRegressor()
In [135]:
model_8.fit(exog_tr_ml, y_train_ml)
Out[135]:
RandomForestRegressor()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
RandomForestRegressor()
In [136]:
param_grid = {
    'n_estimators': [10, 50, 100],
    'max_depth': [10, 20, 30],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4]
}
In [137]:
grid_search = GridSearchCV(estimator=model_8, param_grid=param_grid, cv=5, scoring='neg_mean_squared_error', verbose=2, n_jobs=-1)
In [138]:
best_model = grid_search.fit(exog_tr_ml, y_train_ml)

Fitting 5 folds for each of 81 candidates, totalling 405 fits
In [139]:
print("Mejores parámetros:", grid_search.best_params_)

Mejores parámetros: {'max_depth': 30, 'min_samples_leaf': 1, 'min_samples_split': 2, 'n_estimators': 50}
In [140]:
best_model = RandomForestRegressor(max_depth=30, min_samples_leaf= 1, min_samples_split= 2, n_estimators=50, random_state=23).fit(exog_tr_ml, y_train_ml)
best_pred = best_model.predict(exog_tst_ml)
In [141]:
plt.figure()
plt.plot(y_train_ml.index, y_train_ml, label='Train')
plt.plot(y_test_ml.index, y_test_ml, label='Test')
plt.plot(y_test_ml.index, best_pred, label='Predictions')
plt.xlabel('Date')
plt.ylabel('Value')
plt.title("Modelo 8 - Random Forest"
          "\n MAPE: " + str(round(mean_absolute_percentage_error(y_test_ml, best_pred) * 100)) + "%" +
          "\n R²: " + str(round(r2_score(y_test_ml, best_pred), 2)))
plt.legend()
plt.show()
No description has been provided for this image

EVALUACIÓN MODELO 8

El modelo 6 tiene métricas muy alejadas de los modelos predictivos anteriores. Este modelo no es en absoluto un buen modelo predictivo para la serie temporal objetivo. En la gráfica se puede apreciar que el modelo detecta los flujos estacionales aunque de manera muy poco precisa.

MODELO 9

In [142]:
param_grid = {
    
    'learning_rate': [0.01, 0.1, 0.2],

    'max_depth': [3, 5, 7],

    'subsample': [0.8, 0.9, 1.0]

}

grid_search = GridSearchCV(XGBRegressor(), param_grid, cv=3)

grid_search.fit(exog_tr_ml, y_train_ml)

best_params = grid_search.best_params_
In [143]:
model_9 = xgb.XGBRegressor(**best_params)

model_9.fit(exog_tr_ml, y_train_ml)
Out[143]:
XGBRegressor(base_score=None, booster=None, callbacks=None,
             colsample_bylevel=None, colsample_bynode=None,
             colsample_bytree=None, device=None, early_stopping_rounds=None,
             enable_categorical=False, eval_metric=None, feature_types=None,
             gamma=None, grow_policy=None, importance_type=None,
             interaction_constraints=None, learning_rate=0.1, max_bin=None,
             max_cat_threshold=None, max_cat_to_onehot=None,
             max_delta_step=None, max_depth=7, max_leaves=None,
             min_child_weight=None, missing=nan, monotone_constraints=None,
             multi_strategy=None, n_estimators=None, n_jobs=None,
             num_parallel_tree=None, random_state=None, ...)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
XGBRegressor(base_score=None, booster=None, callbacks=None,
             colsample_bylevel=None, colsample_bynode=None,
             colsample_bytree=None, device=None, early_stopping_rounds=None,
             enable_categorical=False, eval_metric=None, feature_types=None,
             gamma=None, grow_policy=None, importance_type=None,
             interaction_constraints=None, learning_rate=0.1, max_bin=None,
             max_cat_threshold=None, max_cat_to_onehot=None,
             max_delta_step=None, max_depth=7, max_leaves=None,
             min_child_weight=None, missing=nan, monotone_constraints=None,
             multi_strategy=None, n_estimators=None, n_jobs=None,
             num_parallel_tree=None, random_state=None, ...)
In [144]:
predictions = model_9.predict(exog_tst_ml)
In [145]:
plt.figure()
plt.plot(y_train_ml.index, y_train_ml, label='Train')
plt.plot(y_test_ml.index, y_test_ml, label='Test')
plt.plot(y_test_ml.index, predictions, label='Predictions')
plt.xlabel('Date')
plt.ylabel('Value')
plt.title("Modelo 9 - XGB"
          "\n MAPE: " + str(round(mean_absolute_percentage_error(y_test_ml, predictions) * 100)) + "%" +
          "\n R²: " + str(round(r2_score(y_test_ml, predictions), 2)))
plt.legend()
plt.show()
No description has been provided for this image

EVALUACIÓN MODELO 9

El modelo 9 tiene métricas aún más alejadas del modelo anterior. Tampoco este modelo es un buen modelo predictivo para la serie temporal objetivo. En la gráfica se puede apreciar que el modelo detecta los flujos estacionales, aunque de manera muy poco precisa, e ignora del todo la tendencia.

CONCLUSIÓN¶

Con este estudio, en primer lugar, se ha demostrado la importancia de las variable exógenas en la predicción de la serie temporal, puesto que el modelo ARIMA que no cuenta con estas variables (modelo 7), devuelve unas métricas muy escasas.

Cabe añadir que la serie temporal cuenta con más observaciones para los últimos años, con lo cual es probable que este factor haya afectado la capacidad predictiva de algunos modelos.

Los modelos 1, 2 y 5 son los modelos con mejores métricas entre los estudiados, aunque destaca el modelo 5 por su buen coeficiente de determinación y porcentaje de error contenido.

Aquí se detalla una explicación de cada elemento del modelo y de su influencia en el mismo:

  • Coeficientes de las variables exógenas:

    • "Establishments", "Employment", "Unemployment Rate", "Inflation %", "GDP_trillion". Estos coeficientes indican la relación entre estas variables y la variable objetivo en el modelo. Por ejemplo, el coeficiente para "Establishments" es 0.0501, lo que sugiere que un aumento en el número de establecimientos debería estar asociado con un aumento en la variable objetivo, aunque en este caso, los coeficientes de estas variables, no parecen ser significativos, con lo cual no hay evidencia para rechazar la H0, puesto que:
      H0: La variable no tiene efecto significativo

    • "LMA_Central New York Region". Este coeficiente indica la relación con esta variable y la variable objetivo. Una mayor presencia de esta variable sugiere una diminución de la variable objetivo. Esta variable es estadisticamente significativa.

  • Coeficientes MA (Moving Average) (ma.L1, ma.S.L4). Estos coeficientes representan la relación lineal entre los errores pasados y el error actual. En el modelo, el coeficiente del primer retardo regular es estadisticamente significativo, con lo cual se puede rechazar la H0 precedentemente definida.

BIBLIOGRAFÍA

  • https://www.kaggle.com/datasets/thedevastator/new-york-quarterly-employment-and-wage-data-2000/data
  • https://ycharts.com/indicators/us_pce_quarterly_inflation_rate
  • https://data.world/data-ny-gov/5hyu-bdh8
  • https://usafacts.org/metrics/gross-domestic-product-gdp-by-state-new-york/?adjustment=None&timeGranularity=Quarterly
  • https://dol.ny.gov/new-york-state-geography
  • https://dol.ny.gov/lwdbs